Martino Ferrari
1d399ccc9b
retrieve() now takes a CameraModel directly instead of the removed pixel_scale/viewing_angle_deg override kwargs, matching the rest of the pipeline's shared-camera convention. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
105 lines
3.3 KiB
Python
105 lines
3.3 KiB
Python
import numpy as np
|
|
import pytest
|
|
|
|
from he11lib.geometry import CameraModel
|
|
from he11lib.modes import LGBasis
|
|
from he11lib.phase_retrieval import PhaseRetriever, propagate_angular_spectrum
|
|
from he11lib.synthetic import SyntheticBeamGenerator
|
|
|
|
W0 = 5e-3
|
|
Z0 = 0.5
|
|
WAVELENGTH = 1.76e-3
|
|
PIXEL_SCALE = 3e-4
|
|
CAMERA_DISTANCE = 5.0
|
|
IMAGE_SHAPE = (121, 121)
|
|
|
|
|
|
def make_basis():
|
|
return LGBasis(w0=W0, z0=Z0, wavelength=WAVELENGTH)
|
|
|
|
|
|
def make_camera():
|
|
focal_length_px = (CAMERA_DISTANCE + Z0) / PIXEL_SCALE
|
|
return CameraModel(
|
|
focal_length_px=focal_length_px,
|
|
position=(0.0, 0.0, -CAMERA_DISTANCE),
|
|
orientation_deg=(0.0, 0.0, 0.0),
|
|
)
|
|
|
|
|
|
def make_grid():
|
|
coords = (np.arange(IMAGE_SHAPE[0]) - IMAGE_SHAPE[0] // 2) * PIXEL_SCALE
|
|
x, y = np.meshgrid(coords, coords)
|
|
return x, y
|
|
|
|
|
|
def test_propagate_round_trip_recovers_original_field():
|
|
basis = make_basis()
|
|
x, y = make_grid()
|
|
field = basis.field(x, y, Z0, p=0, l=0)
|
|
|
|
forward = propagate_angular_spectrum(field, PIXEL_SCALE, dz=0.05, wavelength=WAVELENGTH)
|
|
back = propagate_angular_spectrum(forward, PIXEL_SCALE, dz=-0.05, wavelength=WAVELENGTH)
|
|
|
|
np.testing.assert_allclose(back, field, atol=1e-3 * np.max(np.abs(field)))
|
|
|
|
|
|
def test_propagate_matches_lgbasis_analytic_evolution():
|
|
basis = make_basis()
|
|
x, y = make_grid()
|
|
field_at_waist = basis.field(x, y, Z0, p=0, l=0)
|
|
|
|
dz = 0.05
|
|
propagated = propagate_angular_spectrum(field_at_waist, PIXEL_SCALE, dz=dz, wavelength=WAVELENGTH)
|
|
analytic = basis.field(x, y, Z0 + dz, p=0, l=0)
|
|
|
|
np.testing.assert_allclose(
|
|
np.abs(propagated) ** 2, np.abs(analytic) ** 2, atol=1e-2 * np.max(np.abs(analytic) ** 2)
|
|
)
|
|
|
|
|
|
def test_retrieve_recovers_pure_mode_purity():
|
|
# Keep z distances close to the waist so the (widening) beam stays well
|
|
# within the frame -- otherwise FFT wraparound/clipping at the edges
|
|
# degrades angular-spectrum propagation accuracy.
|
|
basis = make_basis()
|
|
camera = make_camera()
|
|
gen = SyntheticBeamGenerator(basis=basis, camera=camera)
|
|
z_list = [0.47, 0.5, 0.53]
|
|
planes = gen.generate(
|
|
coefficients={(0, 0): 1.0 + 0j}, z_list=z_list, image_shape=IMAGE_SHAPE, noise_std=1e-5, seed=0
|
|
)
|
|
|
|
retriever = PhaseRetriever(wavelength=WAVELENGTH)
|
|
result = retriever.retrieve(planes, camera, max_iterations=100)
|
|
|
|
dx = float(result.x[0, 1] - result.x[0, 0])
|
|
coeffs = basis.project(
|
|
result.field, result.x, result.y, dx, result.z, modes=[(0, 0), (1, 0), (0, 1)]
|
|
)
|
|
total_power = sum(abs(c) ** 2 for c in coeffs.values())
|
|
purity_00 = abs(coeffs[(0, 0)]) ** 2 / total_power
|
|
assert purity_00 > 0.9
|
|
|
|
|
|
def test_retrieve_estimates_beam_center():
|
|
basis = make_basis()
|
|
camera = make_camera()
|
|
gen = SyntheticBeamGenerator(basis=basis, camera=camera)
|
|
z_list = [0.47, 0.5, 0.53]
|
|
true_center = (15 * PIXEL_SCALE, -8 * PIXEL_SCALE)
|
|
planes = gen.generate(
|
|
coefficients={(0, 0): 1.0 + 0j},
|
|
z_list=z_list,
|
|
image_shape=IMAGE_SHAPE,
|
|
center=true_center,
|
|
noise_std=1e-5,
|
|
seed=1,
|
|
)
|
|
|
|
retriever = PhaseRetriever(wavelength=WAVELENGTH)
|
|
result = retriever.retrieve(planes, camera, max_iterations=100)
|
|
|
|
assert result.center[0] == pytest.approx(true_center[0], abs=3 * PIXEL_SCALE)
|
|
assert result.center[1] == pytest.approx(true_center[1], abs=3 * PIXEL_SCALE)
|