Protoplanetary disk models
============================

:class:`~pyralysis.models.sky.DiskSource` provides **parametric protoplanetary-disk sky
models** for lightweight interferometric workflows. Each component specifies a radial
brightness profile, outer (and optionally inner) radius, inclination, and position angle.
Visibilities are evaluated **directly in the UV plane** using closed-form Fourier transforms
(or a small radial quadrature for power-law disks) — no hydrodynamics and no radiative
transfer.

**What this is for:** fast sky simulation, imaging pipeline tests, AI/ML training data,
and algorithm benchmarking. Compared to hydrodynamic or Monte Carlo radiative-transfer
models, these analytic disks are **lighter, easier to configure, and much faster** to
evaluate (especially at scale with Dask). You trade physical realism for speed and
reproducibility when the goal is to test interferometry software, not to model gas
dynamics.

**What this is not:** a substitute for realistic disk physics (gaps from planet formation,
chemistry, time evolution, etc.). Use parametric rings and envelopes to *mimic* observed
morphology, then validate imaging or ML methods against known ground truth in the
visibility domain. Compose multiple components with
:class:`~pyralysis.models.sky.CompositeSource` (or the ``+`` operator) — see
``examples/notebooks/protoplanetary_disk_alma_simulation.ipynb``.

.. seealso::
   - :doc:`simulation` — general simulation workflow and sky-model overview
   - :doc:`examples` — notebook index including the ALMA disk tutorial
   - :class:`~pyralysis.models.sky.DiskSource` — API reference (auto-generated)

Geometry: inclination and position angle
-----------------------------------------

Disk radii are defined in the **disk midplane**. Inclination :math:`i` and position
angle :math:`\mathrm{PA}` project the disk onto the sky and rotate the UV coordinates
used by the Fourier transform. The formulas below follow
:func:`~pyralysis.models.sky.disk_profiles.projected_uv_radius`.

**Position angle convention:** :math:`\mathrm{PA}` is the major-axis orientation
**east of north** — the same definition as FITS ``BPA``, CASA restoring beams, and
:class:`radio_beam.Beam` (counterclockwise from north on the sky). This matches
ALMA protoplanetary-disk papers. Internally the rotation uses
:math:`\theta = \pi/2 - \mathrm{PA}` because :math:`l` is defined positive toward
east and :math:`m` toward north.

Sky direction-cosine offsets :math:`(l, m)` are in radians (small-angle limit:
:math:`\Delta\mathrm{RA} \approx l / \cos\delta`). :math:`l` is positive toward
**east**; :math:`m` is positive toward **north**. Pyralysis builds image-plane
coordinates like :meth:`~pyralysis.grids.Grid.create_image_meshgrid`:

.. math::

   l = (i - i_0)\,\mathrm{cellsize}_l,\quad
   m = (j - j_0)\,\mathrm{cellsize}_m

with :math:`\mathrm{cellsize}_l = \texttt{cellsize[0]} < 0` and
:math:`\mathrm{cellsize}_m = \texttt{cellsize[1]} > 0`. The negative horizontal
cellsize is what makes :math:`\Delta\mathrm{RA}` run **left to right from
positive to negative** on displayed images. With :math:`\theta = \pi/2 - \mathrm{PA}`,
after rotation into the disk frame,

.. math::

   l_r = l\cos\theta + m\sin\theta,\quad
   m_r = -l\sin\theta + m\cos\theta

.. math::

   \rho = \sqrt{l_r^2 + \left(\frac{m_r}{\cos i}\right)^2}

where :math:`l_r` lies along the projected major axis and :math:`m_r` along the minor
axis. The same :math:`\rho` projection enters the **UV-plane** visibility factor via
:func:`~pyralysis.models.sky.disk_profiles.projected_uv_radius`.

In the UV plane, coordinates :math:`(u, v)` in wavelengths use the same
:math:`\theta` and :math:`\cos i` projection as in the class docstring of
:class:`~pyralysis.models.sky.DiskSource`.

:class:`~pyralysis.models.sky.GaussianSource` uses the same :math:`\mathrm{PA}`
convention for ``theta_angle``.

.. note::
   Bright **ring lobes** in inclined disks trace the **minor axis** (limb brightening),
   not the geometric major axis. When comparing dirty images to model
   :math:`\mathrm{PA}`, measure the **envelope** orientation or use visibility-domain
   geometry checks rather than peak-to-peak separation alone.

Model images (sky brightness)
-----------------------------

Before convolution with the instrumental response, the disk appears on the sky as a
2D brightness distribution :math:`I(l, m)`. The figures below are
**normalized model images** (not dirty maps or simulated visibilities) rasterized
with the Pyralysis ``cellsize`` sign convention above.

.. figure:: _static/disk_models/disk_profiles_faceon.png
   :alt: Face-on disk profile model images
   :width: 100%

   **Face-on** normalized model images for the five built-in profiles
   (:math:`i = 0^\circ`). :math:`\Delta\mathrm{RA}` increases to the left.

.. figure:: _static/disk_models/disk_profiles_inclined.png
   :alt: Inclined disk profile model images
   :width: 100%

   Same profiles at :math:`i = 60^\circ`, :math:`\mathrm{PA} = 30^\circ`. The projected
   ellipse flattens by :math:`\cos i` along the minor axis.

Sky brightness shapes (normalized to peak 1):

.. list-table::
   :header-rows: 1
   :widths: 16 64

   * - Profile
     - :math:`I(\rho)` in the disk midplane (:math:`\rho` in rad)
   * - ``uniform``
     - Constant inside :math:`\rho \le r_\mathrm{out}`; zero outside
   * - ``ring``
     - Constant for :math:`r_\mathrm{in} \le \rho \le r_\mathrm{out}`
   * - ``thin_ring``
     - Infinitely thin shell at :math:`\rho = r_\mathrm{out}`; the documentation
       figure uses a narrow annulus for display only
   * - ``exponential``
     - :math:`\exp(-\rho / r_\mathrm{out})`
   * - ``power_law``
     - :math:`\rho^{-p}` between :math:`r_\mathrm{in}` and :math:`r_\mathrm{out}`
       (peak-normalized for display)

Regenerate the documentation figures after changing profile math (one-off script,
not part of the installed package):

.. code-block:: bash

   python docs/scripts/generate_disk_profile_figures.py

UV-plane visibilities (Fourier theory)
--------------------------------------

:class:`~pyralysis.models.sky.DiskSource` returns visibilities

.. math::

   V(u,v) = I_0 \times F(u,v) \times e^{-2\pi i (u l_0 + v m_0 + w (n_0 - 1))}

where :math:`I_0` is the Stokes flux from ``reference_intensity``,
:math:`F` is the **normalized** spatial Fourier factor from
:func:`~pyralysis.models.sky.disk_profiles.disk_profile_visibility`, and the
phase term shifts the disk when ``sky_position`` is offset from the phase center.

With :math:`x = 2\pi r_\mathrm{out}\,\rho` and :math:`\rho` the projected UV radius
above, the analytic factors implemented in :mod:`pyralysis.models.sky.disk_profiles`
are:

**Uniform disk**

.. math::

   F_\mathrm{uniform}(\rho) = \frac{2 J_1(x)}{x}

**Annular ring** (:math:`r_\mathrm{in} < r_\mathrm{out}`)

.. math::

   F_\mathrm{ring}(\rho) =
   \frac{r_\mathrm{out}^2}{r_\mathrm{out}^2 - r_\mathrm{in}^2}
   \frac{2 J_1(2\pi r_\mathrm{out}\rho)}{2\pi r_\mathrm{out}\rho}
   - \frac{r_\mathrm{in}^2}{r_\mathrm{out}^2 - r_\mathrm{in}^2}
   \frac{2 J_1(2\pi r_\mathrm{in}\rho)}{2\pi r_\mathrm{in}\rho}

**Thin ring** at :math:`r_\mathrm{out}`

.. math::

   F_\mathrm{thin}(\rho) = J_0(2\pi r_\mathrm{out}\rho)

**Exponential disk**

.. math::

   F_\mathrm{exp}(\rho) = \left(1 + x^2\right)^{-3/2}

.. note::
   The disk ``exponential`` profile is **not** :class:`~pyralysis.models.sky.GaussianSource`.
   A Gaussian has image brightness :math:`\propto \exp(-r^2/2\sigma^2)` and UV visibility
   :math:`\exp(-2\pi^2\sigma^2\rho^2)`; the disk exponential has
   :math:`\exp(-\rho/r_\mathrm{out})` in the midplane and
   :math:`F_\mathrm{exp}=(1+x^2)^{-3/2}` above. Use :class:`~pyralysis.models.sky.GaussianSource`
   for elliptical Gaussians; use ``profile_type="exponential"`` for a radial disk envelope.

**Power law** :math:`\rho^{-p}` between :math:`r_\mathrm{in}` and :math:`r_\mathrm{out}`

.. math::

   F_\mathrm{pl}(\rho) = \int_{r_\mathrm{in}}^{r_\mathrm{out}}
   r^{1-p}\, J_0(2\pi r \rho)\, dr \Big/ \int_{r_\mathrm{in}}^{r_\mathrm{out}} r^{1-p}\, dr

The implementation uses a fixed 100-point radial quadrature (see
``_POWER_LAW_QUADRATURE_STEPS`` in :mod:`pyralysis.models.sky.disk_profiles`).

All factors are evaluated in **float32** with Dask block-wise SciPy Bessel calls for
HPC-friendly simulation.

Building individual disk components
-----------------------------------

Angle attributes accept ``Quantity``, ``Angle``, strings, or numeric literals
(radii in **radians**, orientation angles in **degrees** when passed as bare numbers),
matching :class:`~pyralysis.models.sky.GaussianSource`.

**Uniform disk**

.. code-block:: python

   from astropy import units as u
   from pyralysis.models.sky import DiskSource

   inner_disk = DiskSource(
       reference_intensity=0.5 * u.Jy,
       reference_frequency=230e9 * u.Hz,
       sky_position="16h00m00s -30d00m00s",
       profile_type="uniform",
       major_radius=0.12 * u.arcsec,
       inclination=45 * u.deg,
       position_angle=30 * u.deg,
   )

**Bright ring**

.. code-block:: python

   bright_ring = DiskSource(
       reference_intensity=2.5 * u.Jy,
       reference_frequency=230e9 * u.Hz,
       sky_position="16h00m00s -30d00m00s",
       profile_type="ring",
       inner_radius=0.35 * u.arcsec,
       major_radius=0.55 * u.arcsec,
       inclination=45 * u.deg,
       position_angle=30 * u.deg,
   )

**Thin ring** (delta shell in the UV plane)

.. code-block:: python

   thin_ring = DiskSource(
       reference_intensity=1.0 * u.Jy,
       profile_type="thin_ring",
       major_radius=0.5 * u.arcsec,
       inclination=60 * u.deg,
       position_angle=30 * u.deg,
   )

**Exponential envelope**

.. code-block:: python

   envelope = DiskSource(
       reference_intensity=1.0 * u.Jy,
       profile_type="exponential",
       major_radius=0.7 * u.arcsec,
       inclination=45 * u.deg,
       position_angle=30 * u.deg,
   )

**Power-law disk**

.. code-block:: python

   power_disk = DiskSource(
       reference_intensity=1.0 * u.Jy,
       profile_type="power_law",
       inner_radius=0.2 * u.arcsec,
       major_radius=1.0 * u.arcsec,
       power_index=1.5,
       inclination=45 * u.deg,
       position_angle=30 * u.deg,
   )

Composing multi-component disks
-------------------------------

Add components with :class:`~pyralysis.models.sky.CompositeSource` or ``+``.
Negative ``reference_intensity`` on a ring can carve a **parametric gap** without a
dedicated gap profile.

.. figure:: _static/disk_models/disk_composite_example.png
   :alt: Composite parametric disk model image
   :width: 55%

   Example composite: exponential envelope + bright ring − darker ring + compact
   uniform core (same geometry as the ALMA notebook).

.. code-block:: python

   from pyralysis.models.sky import CompositeSource, DiskSource

   common = dict(
       sky_position="16h00m00s -30d00m00s",
       reference_frequency=230e9,
       inclination=45,
       position_angle=30,
   )

   disk_model = CompositeSource(
       sources=[
           DiskSource(
               reference_intensity=1.0,
               profile_type="exponential",
               major_radius="0.7 arcsec",
               **common,
           ),
           DiskSource(
               reference_intensity=2.5,
               profile_type="ring",
               inner_radius="0.35 arcsec",
               major_radius="0.55 arcsec",
               **common,
           ),
           DiskSource(
               reference_intensity=-0.8,
               profile_type="ring",
               inner_radius="0.22 arcsec",
               major_radius="0.32 arcsec",
               **common,
           ),
           DiskSource(
               reference_intensity=0.5,
               profile_type="uniform",
               major_radius="0.12 arcsec",
               **common,
           ),
       ],
       operation="add",
   )

Simulating visibilities
-----------------------

Once the sky model is defined, pass it to :class:`~pyralysis.simulation.Simulator` as
for any other source. Pyralysis evaluates **visibilities analytically in the UV plane**.
The model-image figures earlier show **true sky brightness**; a dirty map is that
brightness convolved with the synthesised beam and affected by UV-coverage sidelobes.

The ALMA notebook (``protoplanetary_disk_alma_simulation.ipynb``) uses the same composite
``disk_model`` geometry as above. Its observation setup is:

- **Array:** ``alma.cycle9.7.cfg`` (43 antennas, long-baseline Cycle 9)
- **Band 6:** 10 channels of 100 MHz from 229.55 to 230.45 GHz (centred on 230 GHz)
- **Time:** 10 minutes total, 2-minute integrations (five time steps)
- **Imaging:** ``cellsize = \mathrm{theo\_resolution}/7``, fixed ``imsize = 2048``
  (use 1024 for a faster run), Natural weighting, then :class:`~pyralysis.transformers.DirtyMapper`
- **Display:** crop ±1 arcsec around the phase centre so rings are easy to see

Omit explicit ``chunks`` on :class:`~pyralysis.simulation.Simulator` so Dask uses
``"auto"``. Expect similar morphology between model images and the dirty map, not a
pixel-identical match.

.. tip::
   Resolved structure needs radii comparable to or larger than the synthesized beam.
   With ALMA Band 6 and long-baseline configurations the beam can be
   :math:`\sim 0.07''` or finer; use arcsecond-scale radii in tutorials unless you
   deliberately want an unresolved component.

.. code-block:: python

   from astropy import units as u
   from astropy.coordinates import Angle
   from astropy.time import Time

   from pyralysis.io.antenna_config_io import AntennaConfigurationIo
   from pyralysis.simulation import Simulator

   interferometer = AntennaConfigurationIo(
       input_name="src/pyralysis/simulation/antenna_configs/alma.cycle9.7.cfg"
   ).read()
   interferometer.configure_observation(
       min_frequency=229.55e9,
       max_frequency=230.55e9,  # half-open grid → 10 × 100 MHz channels to 230.45 GHz
       frequency_step=100e6,
       right_ascension=Angle("16h00m00s", unit=u.hourangle),
       declination=Angle("-30d00m00s"),
       integration_time=2 * u.minute,
       observation_time="10m",
       reference_time=Time("2024-01-01T00:00:00", format="isot"),
   )

   dataset = Simulator(interferometer=interferometer, sources=disk_model).simulate(
       create_dataset=True
   )

The notebook continues with :class:`~pyralysis.injectors.ThermalNoiseInjector`
(``channel_bandwidth=100 MHz``, ``integration_time=2 min``), Natural weighting, and
dirty imaging. See :doc:`simulation` for noise injection, chunking, and dataset export.

See also
--------

- :doc:`simulation`
- :doc:`examples`
- :doc:`faraday` — Faraday rotation for polarized components

----

:doc:`simulation` | :doc:`examples`