.. image:: https://mybinder.org/badge_logo.svg
:target: https://mybinder.org/v2/gh/ehpor/hcipy/HEAD?labpath=doc/tutorial_notebooks/AtmosphericSimulation/AtmosphericSimulation.ipynb
:alt: Open in Binder
Simulating Atmospheric Turbulence
=================================
We’ll model how Earth’s atmosphere distorts starlight using HCIPy’s
phase screen framework, from single-layer models to realistic
multi-layer atmospheric profiles.
Atmospheric turbulence causes spatial and temporal variations in the
refractive index along the line of sight, which translates to optical
path length differences — phase aberrations — that degrade image
quality. In this tutorial we’ll build up from a single turbulent layer
to a full multi-layer profile, exploring how parameters like the Fried
parameter :math:`r_0` and wind speed affect the Point Spread Function
(PSF).
.. code:: ipython3
import hcipy
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
We start by defining the telescope and observing parameters. The
aperture diameter sets the diffraction limit, while the wavelength
determines the spatial scale of the focal plane. We’ll use a
:math:`256 \times 256` pupil grid with generous padding to avoid
aliasing, and a focal grid that oversamples the PSF by a factor of 4
(:math:`q = 4`) to resolve fine structure.
.. code:: ipython3
# Telescope parameters
pupil_diameter = 8.2 # meters
wavelength = 1.5e-6 # meters
# Computational grids
pupil_grid = hcipy.make_pupil_grid(256, pupil_diameter * 1.2)
spatial_resolution = wavelength / pupil_diameter
focal_grid = hcipy.make_focal_grid(q=4, num_airy=32, spatial_resolution=wavelength / pupil_diameter)
Now that we have our grids defined, we can set up the Fraunhofer
propagator and the telescope aperture. The propagator will transform the
electric field from the pupil plane to the focal plane.
.. code:: ipython3
# Propagator and aperture
prop = hcipy.FraunhoferPropagator(pupil_grid, focal_grid)
aperture = hcipy.make_vlt_aperture()(pupil_grid)
print(f"Telescope diameter: {pupil_diameter} m")
print(f"Observing wavelength: {wavelength*1e9:.0f} nm")
print(f"Pupil grid: {pupil_grid.dims} points")
print(f"Focal grid: {focal_grid.dims} points")
.. parsed-literal::
Telescope diameter: 8.2 m
Observing wavelength: 1500 nm
Pupil grid: (256, 256) points
Focal grid: (256, 256) points
The core of an atmospheric simulation is the phase screen — a 2D map of
the optical path difference introduced by refractive index fluctuations.
We’ll use HCIPy’s ``InfiniteAtmosphericLayer``, which generates
scrollable turbulence patterns using the method of Assemat et
al. (2006). This type of layer synthesizes an infinitely large phase
screen on-the-fly by summing Fourier modes with Kolmogorov-von Kármán
statistics, so wind can blow the screen across the aperture without
repeating.
The key parameters are the Fried parameter :math:`r_0` (which sets the
turbulence strength), the outer scale :math:`L_0` (the largest eddy
size), and the wind velocity (which controls the frozen-flow evolution
rate).
.. code:: ipython3
# Atmospheric parameters
fried_parameter = 0.2 # r0 = 20 cm
outer_scale = 20 # L0 (meters)
velocity = 10 # wind speed (m/s)
# Convert Fried parameter to Cn squared
Cn_squared = hcipy.Cn_squared_from_fried_parameter(fried_parameter, wavelength=500e-9)
print("Atmospheric Parameters:")
print(f" Fried parameter (r0): {fried_parameter} m at 500 nm")
print(f" Outer scale (L0): {outer_scale} m")
print(f" Wind velocity: {velocity} m/s")
print(f" Cn squared: {Cn_squared:.2e} m^(-2/3)")
.. parsed-literal::
Atmospheric Parameters:
Fried parameter (r0): 0.2 m at 500 nm
Outer scale (L0): 20 m
Wind velocity: 10 m/s
Cn squared: 2.19e-13 m^(-2/3)
With these parameters defined, we can create the
``InfiniteAtmosphericLayer``. The layer will use these atmospheric
parameters to generate phase screens consistent with Kolmogorov-von
Kármán turbulence statistics.
.. code:: ipython3
# Atmospheric layer
layer = hcipy.InfiniteAtmosphericLayer(
pupil_grid,
Cn_squared=Cn_squared,
L0=outer_scale,
velocity=velocity,
seed=0
)
print("\nAtmospheric phase screen created")
.. parsed-literal::
Atmospheric phase screen created
Let’s take an instantaneous snapshot of the phase screen and compute
wavefront error statistics across the aperture.
.. code:: ipython3
# Phase screen at observing wavelength
phase = layer.phase_for(wavelength)
# Visualize
plt.figure(figsize=(10, 8))
im = hcipy.imshow_field(phase, cmap='RdBu_r', mask=aperture, vmin=-10, vmax=10)
plt.colorbar(im, label='Phase (radians)')
plt.title(f'Atmospheric Phase Screen\nr0 = {fried_parameter}m at {wavelength*1e9:.0f}nm')
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.axis('equal')
plt.tight_layout()
plt.show()
# Wavefront error statistics
phase_inside_aperture = phase[aperture > 0]
rms_wfe = np.sqrt(np.mean(phase_inside_aperture**2))
peak_to_valley = np.max(phase_inside_aperture) - np.min(phase_inside_aperture)
print(f"Wavefront Error Statistics (inside aperture):")
print(f" RMS: {rms_wfe:.2f} rad ({rms_wfe * wavelength * 1e9:.1f} nm)")
print(f" Peak-to-Valley: {peak_to_valley:.2f} rad ({peak_to_valley * wavelength * 1e9:.1f} nm)")
.. image:: output_11_0.png
.. parsed-literal::
Wavefront Error Statistics (inside aperture):
RMS: 2.83 rad (4247.1 nm)
Peak-to-Valley: 12.95 rad (19418.7 nm)
The RMS wavefront error is on the order of several radians —
significantly larger than the :math:`\lambda/14` (≈0.45 rad) Marechal
criterion for diffraction-limited imaging. These phase errors scatter
light from the core of the PSF into a broad halo. Let’s quantify this by
comparing the diffraction-limited PSF to the turbulence-degraded
version.
.. code:: ipython3
# Plane wave entering the telescope
wf = hcipy.Wavefront(aperture, wavelength)
# Apply atmosphere and propagate
wf_atmos = layer(wf)
img_turbulent = prop(wf_atmos)
# Diffraction-limited comparison
img_diffraction_limited = prop(wf)
# Visualize
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Diffraction-limited PSF
im0 = hcipy.imshow_field(
np.log10(img_diffraction_limited.intensity / img_diffraction_limited.intensity.max()),
vmin=-5, vmax=0, cmap='inferno', grid_units=spatial_resolution, ax=axes[0]
)
axes[0].set_title('Diffraction-Limited (No Atmosphere)', fontsize=12)
axes[0].set_xlabel('x (lambda/D)')
axes[0].set_ylabel('y (lambda/D)')
plt.colorbar(im0, ax=axes[0], label='log10(Intensity)')
# Turbulent PSF
im1 = hcipy.imshow_field(
np.log10(img_turbulent.intensity / img_turbulent.intensity.max()),
vmin=-5, vmax=0, cmap='inferno', grid_units=spatial_resolution, ax=axes[1]
)
axes[1].set_title(f'With Atmospheric Turbulence (r0 = {fried_parameter}m)', fontsize=12)
axes[1].set_xlabel('x (lambda/D)')
axes[1].set_ylabel('y (lambda/D)')
plt.colorbar(im1, ax=axes[1], label='log10(Intensity)')
plt.tight_layout()
plt.show()
# Strehl ratio
strehl_dl = np.max(img_diffraction_limited.intensity) / np.sum(img_diffraction_limited.intensity)
strehl_turb = np.max(img_turbulent.intensity) / np.sum(img_turbulent.intensity)
print(f"\nImage Quality Comparison:")
print(f" Diffraction-limited peak intensity: {strehl_dl:.4f}")
print(f" Turbulent peak intensity: {strehl_turb:.4f}")
print(f" Strehl ratio: {strehl_turb/strehl_dl:.2%}")
.. image:: output_13_0.png
.. parsed-literal::
Image Quality Comparison:
Diffraction-limited peak intensity: 0.0460
Turbulent peak intensity: 0.0034
Strehl ratio: 7.41%
Even with :math:`r_0 = 20` cm — which is moderate seeing — the Strehl
ratio drops to just a few percent. This makes it clear why every major
observatory needs adaptive optics to recover diffraction-limited
performance.
The atmosphere isn’t static. Wind moves turbulent layers across the
aperture, causing the PSF to evolve on millisecond timescales. We’ll
simulate this frozen-flow evolution by advancing the
``InfiniteAtmosphericLayer`` in time and recording frames.
.. code:: ipython3
# Animation parameters
duration = 1.0 # seconds
fps = 40
num_frames = int(duration * fps)
time_step = duration / num_frames
writer = hcipy.FFMpegWriter('atmosphere_evolution.mp4', framerate=10)
print(f"Creating animation with {num_frames} frames...")
.. parsed-literal::
Creating animation with 40 frames...
Now we’ll loop through the time steps, evolving the atmosphere,
computing the PSF, and recording frames to a video file. At each step we
display the phase screen alongside the PSF with the current Strehl ratio
overlaid.
.. code:: ipython3
for t in tqdm(np.arange(num_frames) * time_step):
layer.evolve_until(t)
wf_atmos = layer(wf)
img = prop(wf_atmos)
phase = layer.phase_for(wavelength)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Phase pattern
im0 = hcipy.imshow_field(phase, cmap='RdBu_r', mask=aperture, vmin=-10, vmax=10, ax=axes[0])
axes[0].set_title(f'Phase Screen (t = {t*1000:.0f} ms)', fontsize=12)
axes[0].set_xlabel('x (m)')
axes[0].set_ylabel('y (m)')
plt.colorbar(im0, ax=axes[0], label='Phase (rad)')
# Point spread function
im1 = hcipy.imshow_field(np.log10(img.intensity / img.intensity.max()), vmin=-4, cmap='inferno', grid_units=spatial_resolution, ax=axes[1])
axes[1].set_title('Point Spread Function', fontsize=12)
axes[1].set_xlabel('x (lambda/D)')
axes[1].set_ylabel('y (lambda/D)')
plt.colorbar(im1, ax=axes[1], label='log10(Intensity)')
# Strehl annotation
strehl = np.max(img.intensity) / np.sum(img.intensity)
fig.suptitle(f'Strehl Ratio: {strehl/strehl_dl:.1%}', fontsize=14, y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.96])
writer.add_frame(fig)
plt.close(fig)
writer.close()
writer
.. parsed-literal::
100%|██████████| 40/40 [00:32<00:00, 1.25it/s]
.. raw:: html
The animation shows the phase screen drifting across the pupil while the
PSF core jumps and speckles shift. Notice how the Strehl ratio
fluctuates frame-to-frame — it can vary by a factor of 2 or more as
different turbulent eddies cross the aperture.
In reality, atmospheric turbulence is concentrated in discrete layers at
different altitudes — the boundary layer, the jet stream, and the
tropopause — each with its own wind speed and direction. HCIPy provides
``make_mauna_kea_atmospheric_layers`` which returns a realistic profile
based on measurements at Mauna Kea. We’ll randomize the wind directions
of each layer to produce a more realistic, non-degenerate speckle
pattern, then combine them into a ``MultiLayerAtmosphere``.
.. code:: ipython3
# Mauna Kea atmospheric profile
layers = hcipy.make_mauna_kea_atmospheric_layers(pupil_grid, outer_scale=outer_scale)
# Randomize layer orientations
for layer in layers:
angle = np.random.uniform(0, 2 * np.pi)
velocity = np.linalg.norm(layer.velocity)
layer.velocity = [np.cos(angle) * velocity, np.sin(angle) * velocity]
# Multi-layer atmosphere
atmos = hcipy.MultiLayerAtmosphere(layers, scintillation=True)
print(f"Mauna Kea Atmospheric Profile:")
print(f" Number of layers: {len(layers)}")
print(f"\nLayer details:")
for i, layer in enumerate(layers):
print(f" Layer {i+1}: {layer.height:5.0f}m altitude, "
f"Cn squared={layer.Cn_squared:.2e}, velocity={layer.velocity} m/s")
.. parsed-literal::
Mauna Kea Atmospheric Profile:
Number of layers: 6
Layer details:
Layer 1: 500m altitude, Cn squared=5.00e-14, velocity=[ 5.65557002 -8.24709208] m/s
Layer 2: 1000m altitude, Cn squared=1.93e-14, velocity=[-9.35552189 -3.53188478] m/s
Layer 3: 2000m altitude, Cn squared=1.46e-14, velocity=[ 9.84856669 -1.73370535] m/s
Layer 4: 4000m altitude, Cn squared=3.19e-14, velocity=[-6.14066207 -7.89254517] m/s
Layer 5: 8000m altitude, Cn squared=7.34e-14, velocity=[-5.5589905 -8.3124981] m/s
Layer 6: 16000m altitude, Cn squared=2.96e-14, velocity=[ 9.51284607 -3.08314119] m/s
The Mauna Kea profile includes layers spanning altitudes from a few
hundred meters to nearly 20 km, with the strongest turbulence in the
boundary layer and the jet stream. Each layer contributes different
amounts to the total wavefront error.
With a single thin layer, the atmosphere only imprints phase errors. But
when light propagates through multiple layers at different altitudes,
diffraction converts phase variations into amplitude variations — this
is scintillation (the same phenomenon that makes stars twinkle). HCIPy’s
``MultiLayerAtmosphere`` handles this Fresnel propagation between layers
when ``scintillation=True``. The scintillation index — the normalized
variance of the intensity — quantifies how strongly the atmosphere
modulates the amplitude of the wavefront.
.. code:: ipython3
# Set turbulence strength
atmos.Cn_squared = hcipy.Cn_squared_from_fried_parameter(0.1, wavelength=550e-9)
atmos.reset()
# Input wavefront
wf_in = hcipy.Wavefront(pupil_grid.ones(), wavelength=wavelength)
We’ve set up a strong turbulence case (:math:`r_0 = 10` cm) and a
unit-amplitude input wavefront. Now let’s animate the evolution of both
the amplitude (scintillation) and phase patterns over one second.
.. code:: ipython3
# Animation
duration = 1 # seconds
fps = 40
num_frames = int(duration * fps)
time_step = duration / num_frames
writer = hcipy.FFMpegWriter('scintillation_evolution.mp4', framerate=10)
print(f"Creating scintillation animation with {num_frames} frames...")
.. parsed-literal::
Creating scintillation animation with 40 frames...
The animation loop evolves the multi-layer atmosphere, computes the
resulting wavefront, and records the intensity and phase patterns along
with the scintillation index.
.. code:: ipython3
for t in tqdm(np.arange(num_frames) * time_step):
atmos.evolve_until(t)
wf_out = atmos(wf_in)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Intensity (scintillation)
im0 = hcipy.imshow_field(wf_out.intensity * aperture, cmap='gray', ax=axes[0], vmin=0, vmax=3)
axes[0].set_title(f'Intensity (t = {t*1000:.0f} ms)', fontsize=12)
axes[0].set_xlabel('x (m)')
axes[0].set_ylabel('y (m)')
plt.colorbar(im0, ax=axes[0], label='Intensity')
# Phase
im1 = hcipy.imshow_field(wf_out.phase * aperture, cmap='RdBu_r', ax=axes[1], vmin=-3, vmax=3)
axes[1].set_title('Phase', fontsize=12)
axes[1].set_xlabel('x (m)')
axes[1].set_ylabel('y (m)')
plt.colorbar(im1, ax=axes[1], label='Phase (rad)')
# Scintillation statistics
intensity_in_aperture = wf_out.intensity[aperture > 0]
scintillation_index = np.std(intensity_in_aperture) / np.mean(intensity_in_aperture)
fig.suptitle(f'Scintillation Index: {scintillation_index:.3f}', fontsize=14, y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.96])
writer.add_frame(fig)
plt.close(fig)
writer.close()
writer
.. parsed-literal::
100%|██████████| 40/40 [00:25<00:00, 1.56it/s]
.. raw:: html
The intensity across the pupil fluctuates significantly as different
turbulent cells drift through the beam. A scintillation index above 0.1
indicates strong amplitude modulation, which conventional adaptive
optics cannot correct with a deformable mirror alone (it requires a
second DM for amplitude control).
We can also sweep across different seeing conditions without recreating
the layers. The ``atmos.Cn_squared`` property scales the turbulence
strength of all layers simultaneously. Let’s compare three
representative seeing values: 1.0 arcsec (poor), 0.5 arcsec (median),
and 0.25 arcsec (excellent).
.. code:: ipython3
# Reset atmosphere
atmos.reset()
# Seeing conditions
seeing_arcsec = [1.0, 0.5, 0.25]
We’ll create a figure with three subplots and populate each one with the
PSF for a different seeing value, annotating each with the Strehl ratio.
.. code:: ipython3
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
for i, seeing in enumerate(seeing_arcsec):
r0 = hcipy.seeing_to_fried_parameter(seeing)
atmos.Cn_squared = hcipy.Cn_squared_from_fried_parameter(r0)
atmos.reset()
wf_out = atmos(wf_in)
img = prop(wf_out)
im = hcipy.imshow_field(
np.log10(img.intensity / img.intensity.max()),
vmin=-3, cmap='inferno', grid_units=spatial_resolution, ax=axes[i]
)
axes[i].set_title(f'r0 = {r0*100:.0f} cm at 500nm\n({seeing:.2f} arcsec seeing)', fontsize=11)
axes[i].set_xlabel('x (lambda/D)')
axes[i].set_ylabel('y (lambda/D)')
# Strehl ratio
strehl = np.max(img.intensity) / np.sum(img.intensity)
axes[i].text(0.02, 0.98, f'Strehl: {strehl:.2%}',
transform=axes[i].transAxes, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
plt.suptitle('Effect of Seeing on Image Quality', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
.. image:: output_29_0.png
Better seeing (larger :math:`r_0`) produces tighter PSFs with higher
Strehl ratios — at 0.25 arcsec the Strehl is nearly 5× better than at
1.0 arcsec. Poor seeing smears the image into a broad halo with little
distinguishable structure. The multi-layer model also produces more
realistic speckle patterns than a single layer, with the characteristic
broken spoke-like structure seen in real coronagraphic images.
.. code:: ipython3
# Cleanup
import os
if os.path.exists('atmosphere_evolution.mp4'):
os.remove('atmosphere_evolution.mp4')
if os.path.exists('scintillation_evolution.mp4'):
os.remove('scintillation_evolution.mp4')