Optilab

Optilab is a MATLAB library for computational wave optics based on Fourier optics and scalar diffraction theory. It provides:

• a numerical propagation engine (Fresnel, Fraunhofer, angular spectrum),
• physically grounded optical source models (point source, Gaussian beam),
• a Wave representation carrying amplitude, phase, polarization, and sampling,
• a base class for optical elements that modify a field.

The library is written in MATLAB using SI units throughout.

Status: research-grade / in active development. APIs may change.

---

Table of Contents

• Overview

• Features

• Architecture

  • 1. Core Layer (Physics and Numerics)
  • 2. Object Model (Fields, Sources, Elements)
  • 3. Visualization and UI Layer

• Directory Structure

• Example Usage

• Development Roadmap

• Development Guidelines

• License

• Author and Contact

---

Overview

Optilab is designed for simulation of coherent optical wavefields and phenomena such as:

• free-space diffraction,
• propagation of Gaussian beams,
• interference between coherent sources,
• spherical wave emission from point sources.

The numerical propagation models implement scalar solutions to the Helmholtz equation in homogeneous media, with optional paraxial approximations where appropriate.

The long-term goal is to support both:

• script-driven optical experiments (teaching, research, rapid prototyping),
• an interactive optical bench GUI for assembling sources and elements and inspecting fields.

All code and documentation are in English. All physical quantities are SI. Core constants follow CODATA 2018 values [CODATA 2018].

---

Features

Implemented:

• Full propagation engine

  • Automatic selection between Fraunhofer (far field), Fresnel (paraxial), and Angular Spectrum (exact scalar Helmholtz).
  • Forward and backward propagation over arbitrary distances dz.
  • Works in any homogeneous, isotropic medium with refractive index n.

• Physical source models

  • CoherentGaussianBeamSource: TEM(_{00}) Gaussian beam with waist, Gouy phase, curvature, and polarization.
  • CoherentPointSource: ideal isotropic coherent point emitter producing a spherical wave with (1/r) decay and physically scaled amplitude based on radiant power (P) [W].
  • All coherent sources emit a Wave object on an arbitrary observation plane.

• Wave representation

  • Wave: amplitude, phase, polarization (Jones vector), wavelength, spatial sampling, and axial plane z.

  • SphericalWave: concrete subclass of Wave built from a point source model.

  • Convenience methods:

    • getFieldComponents() → complex field components (U_x, U_y)
    • getFieldStack() → Ny × Nx × 2 complex array
    • getIntensity() → irradiance in W/m²
    • getPower() → total optical power via numerical integration

• Polarization support

  • Fields can be stored and propagated either as a scalar complex field or as a 2-component Jones field ([..., :, 1] = Ux, [..., :, 2] = Uy).
  • OpticalElement.apply() can operate either as a scalar transmission mask or as a spatially varying 2×2 Jones matrix.

• Physical constants

  • core/constants.m provides (c), (\varepsilon_0), (\mu_0), etc., in SI units [CODATA 2018].
  • Used for radiometric scaling (e.g. converting source power to electric field amplitude and computing irradiance).

Planned / in progress:

• Thin optical elements such as lenses, apertures, and polarizers as subclasses of OpticalElement.
• Chaining of elements and propagation into a full optical system model.
• High-level plotting utilities and GUI tooling.

---

Architecture

The library is organized in layered fashion to keep physics, object modeling, and visualization separated.

1. Core Layer (Physics and Numerics)

Located in core/.

Responsibilities:

• Propagation of a complex field through space.
• Frequency grid generation and sampling utilities.
• Physical constants.

Key components:

• propagate.m High-level propagation routine. Chooses automatically between:

  • Fraunhofer (far-field / Fourier plane),
  • Fresnel (paraxial),
  • Angular Spectrum (exact scalar Helmholtz, including evanescent components),
  • Identity (dz = 0).

  The selection is based on:

  • propagation distance vs. estimated Fraunhofer distance (z_\mathrm{FF} \approx 2 a^2 / \lambda),
  • numerical aperture content estimated from local phase gradients in the field,
  • basic paraxial validity checks.

  Returns the propagated complex field, the output sampling coordinates, and a diagnostic info struct (method used, NA estimate, etc.).

• propagateFraunhofer.m Explicit Fraunhofer diffraction. Produces the field in the far field / focal plane, and returns the corresponding physical coordinates in the observation plane.

• propagateFresnel.m FFT-based Fresnel propagation using the quadratic phase transfer function in spatial-frequency space.

• propagateAngularSpectrum.m Angular spectrum method. Decomposes the input field into plane waves via a centered 2D FFT, applies (\exp(i k_z dz)) to each spatial frequency component, and reconstructs the propagated field via inverse FFT. This is the most general scalar propagator and is valid beyond the paraxial regime.

• utils/fft2shift.m, utils/ifft2shift.m Centered forward / inverse 2D Fourier transforms.

• utils/makeFrequencyGrid.m Generates spatial frequency axes (in cycles/m) consistent with fftshift ordering, given a sampled spatial grid.

• constants.m Returns a struct of physical constants in SI units (speed of light, vacuum permittivity, impedance of free space, etc.) [CODATA 2018].

Assumptions of the core layer:

• Homogeneous, isotropic, non-absorbing medium of refractive index n (constant over dz).
• Scalar diffraction model (no vector diffraction, no birefringence).
• Periodicity / wrap-around effects can occur if the sampled field is not sufficiently padded before propagation (standard FFT-based limitation).

2. Object Model (Fields, Sources, Elements)

Located in objects/.

This layer provides a physically meaningful interface for constructing optical scenes.

Wave-level classes

• Wave Represents a monochromatic optical field on a transverse plane (z = \text{constant}). Stores:

  • amplitude (A(x,y)) [V/m],
  • phase (\phi(x,y)) [rad],
  • wavelength (\lambda) [m],
  • normalized Jones polarization vector,
  • sampling grid (x, y),
  • plane position z.

  Also provides:

  • irradiance (I(x,y)) in W/m² using ( I = \tfrac{1}{2} c \varepsilon_0 (|U_x|^2 + |U_y|^2) ) [CODATA 2018],
  • total optical power via numerical integration over the plane.

• SphericalWave Subclass of Wave. Builds the field produced by an isotropic coherent point source of radiant power (P) at position ((x_0,y_0,z_0)), evaluated on an arbitrary observation plane (z = z_\mathrm{plane}). Uses ( U(r) = A_0 \exp(i k r) / r ) with ( |A_0| \propto \sqrt{P / (c \varepsilon_0)} ) so that radiometric power is consistent.

Source hierarchy

• Source (abstract) General contract for anything that can emit a field. A source knows:

  • its spatial position [x0 y0 z0] [m],
  • its emission wavelength [m],
  • its polarization (Jones vector),
  • optional metadata (Name, Notes).

  It must implement:

  wave = emit(obj, x, y, z_plane)

  which returns a Wave sampled on the requested grid at axial position z_plane.

• CoherentSource (abstract) Specialization of Source for fully coherent, monochromatic emitters. Adds:

  • deterministic phase,
  • FieldPhaseOffset (global phase offset [rad]) for controlled interference between multiple coherent sources.

  Any subclass of CoherentSource returns a deterministic Wave (not an ensemble average).

• CoherentPointSource Coherent point emitter with total radiant power Power [W]. Generates a SphericalWave on demand. The field decays as (1/r) and accumulates phase (k r).

• CoherentGaussianBeamSource TEM(_{00})-like Gaussian beam under the standard paraxial Gaussian beam model:

  • waist radius (w_0),
  • Rayleigh range (z_R = \pi w_0^2 / \lambda),
  • radius of curvature (R(z)),
  • Gouy phase (\psi(z)),
  • spot size (w(z)),
  • user-defined on-axis amplitude at the waist,
  • uniform Jones polarization.

  emit(...) evaluates the beam on an arbitrary plane z_plane and returns a Wave.

Optical elements

• OpticalElement (abstract) Base class for thin elements that apply a local complex transmission on a plane:

  • scalar transmission mask (T(x,y)), or
  • spatially varying Jones matrix (J(x,y)) acting on ([U_x; U_y]).

  Provides:

  Uout = obj.apply(Uin, x, y, wavelength)

  which applies the element to an incident field (scalar or Jones 2-component).

Planned subclasses include thin lenses, apertures, slits, spatial filters, and polarizing elements. These will inherit from OpticalElement and implement a physically meaningful transmission(...) model (e.g. quadratic phase for an ideal thin lens).

3. Visualization and UI Layer

Located in ui/ (placeholder).

Planned responsibilities:

• Field inspection tools (intensity maps, phase maps, log-scale plots).
• Simple diagnostic inspectors for beams and elements.
• A MATLAB App Designer GUI acting as an “optical bench”: arrange a source, place elements at known z positions, propagate between planes, and visualize the field on detectors.

This layer is not yet finalized.

---

Directory Structure

Optilab/
│
├── core/                          % Physical engine and numerical tools
│   ├── constants.m                % CODATA 2018 physical constants [CODATA 2018]
│   ├── propagation/
│   │   ├── propagate.m            % High-level adaptive propagator
│   │   ├── propagateFresnel.m     % Fresnel (paraxial) propagation
│   │   ├── propagateFraunhofer.m  % Fraunhofer (far-field) propagation
│   │   └── propagateAngularSpectrum.m % Angular spectrum (exact scalar Helmholtz)
│   └── utils/
│       ├── fft2shift.m
│       ├── ifft2shift.m
│       └── makeFrequencyGrid.m
│
├── objects/                       % Physical abstractions
│   ├── Wave.m
│   ├── SphericalWave.m
│   ├── Source.m
│   ├── CoherentSource.m
│   ├── CoherentPointSource.m
│   ├── CoherentGaussianBeamSource.m
│   └── OpticalElement.m           % Base class for thin optical elements
│                                  % (future: Lens, Aperture, Polarizer, etc.)
│
├── ui/                            % Visualization / GUI layer (planned)
│   └── ...                        % App Designer tooling, plot helpers
│
├── examples/                      % Demonstration and validation scripts
│   ├── propagate_test.mlx
│   ├── source_test.mlx
│   ├── Young.mlx
│   └── wave_test.mlx
│
└── README.md

---

Example Usage

Below is a minimal end-to-end example:

1. define a spatial grid,
2. create a coherent Gaussian beam source,
3. emit the field as a Wave at (z = 0),
4. propagate it by 20 cm in air using the adaptive propagator,
5. compute and plot the intensity on the observation plane.

% -------------------------------------------------------------------------
% 1. Spatial sampling grid (meters)
% -------------------------------------------------------------------------
x = linspace(-2e-3, 2e-3, 512);      % 4 mm wide window, 512 samples
y = linspace(-2e-3, 2e-3, 512).';    % column vector for y

% -------------------------------------------------------------------------
% 2. Define a coherent Gaussian beam source
%    Waist at z = 0, linearly polarized along x̂
% -------------------------------------------------------------------------
src = CoherentGaussianBeamSource( ...
    'Waist',         0.5e-3, ...      % beam waist radius w0 [m]
    'PeakAmplitude', 1.0,    ...      % on-axis |E| at waist [V/m]
    'Wavelength',    633e-9, ...      % 633 nm
    'Position',      [0 0 0], ...     % waist located at (0,0,0)
    'Polarization',  [1; 0], ...
    'Name',          "GaussianBeam");

% -------------------------------------------------------------------------
% 3. Emit the field on the z = 0 plane
%    wave0 is a Wave instance
% -------------------------------------------------------------------------
z0 = 0.0;
wave0 = src.emit(x, y, z0);

% Extract complex field components Ux, Uy [V/m]
U0 = wave0.getFieldStack();   % Ny-by-Nx-by-2

% -------------------------------------------------------------------------
% 4. Propagate the field forward by dz = 0.20 m in air (n = 1)
%    The propagate() function will choose Fresnel / Angular Spectrum / etc.
% -------------------------------------------------------------------------
dz = 0.20;     % meters
n_medium = 1;  % air

[U1, x1, y1, info] = propagate( ...
    U0, ...
    wave0.x, ...
    wave0.y, ...
    wave0.Wavelength, ...
    n_medium, ...
    dz);

% U1 is Ny-by-Nx-by-2 complex at z = 0.20 m
% x1, y1 are the sampling coordinates of that output plane.
% info.MethodUsed reports which propagator was selected.

% -------------------------------------------------------------------------
% 5. Compute intensity [W/m^2] and display
%    I = 0.5 * c * eps0 * (|Ux|^2 + |Uy|^2)
% -------------------------------------------------------------------------
const = constants();
Ux1 = U1(:,:,1);
Uy1 = U1(:,:,2);

I = 0.5 * const.c * const.epsilon0 .* (abs(Ux1).^2 + abs(Uy1).^2);

figure;
imagesc(x1*1e3, y1*1e3, I);   % plot in mm
axis image;
xlabel('x [mm]');
ylabel('y [mm]');
title('Intensity at z = 0.20 m');
colorbar;

This example demonstrates:

• source modeling through CoherentGaussianBeamSource,
• physically scaled field generation via Wave,
• propagation using propagate,
• intensity calculation using physical units.

---

Development Roadmap

Phase	Objective	Status / Notes
1	Core numerical propagation engine	Done. Fresnel / Fraunhofer / Angular Spectrum implemented.
2	Coherent source and field model	Done. Wave, SphericalWave, CoherentPointSource, CoherentGaussianBeamSource.
3	Optical elements	In progress. OpticalElement base class exists (scalar mask / Jones matrix). Thin lens, aperture, polarizer, etc. to be added as concrete subclasses.
4	System-level chaining	Planned. High-level orchestration of: emit → apply element(s) → propagate → detect. (An OpticalSystem manager class will live here.)
5	Interference / holography scenarios	Planned. Multi-source interference, Young-type experiments, off-axis holography.
6	Visualization / GUI	Planned. App Designer interface for interactive optical bench, plotting of amplitude/phase/intensity across planes.

---

Development Guidelines

• Language All code, comments, commit messages, and documentation are in English.

• Units and conventions

  • All physical quantities are SI.
  • Positions are in meters.
  • Wavelength is in meters.
  • Electric field amplitude is in V/m.
  • Intensities are in W/m².
  • Propagation assumes the optical axis is +z.

• Sampling conventions

  • x is a row vector (1 × Nx) and corresponds to columns.
  • y is a column vector (Ny × 1) and corresponds to rows.
  • (x,y) define a transverse plane at axial position z.
  • Frequency grids (fx, fy) follow fftshift ordering.

• Polarization

  • Polarization is expressed as a 2-element Jones vector [px; py].
  • Wave stores normalized polarization.
  • A propagated field can be represented as an Ny × Nx × 2 array [Ux, Uy].

• Coding style

  • Classes: UpperCamelCase (e.g. CoherentGaussianBeamSource, OpticalElement).
  • Methods / variables: lowerCamelCase (e.g. emit, getFieldStack, makeFrequencyGrid).
  • Constants: ALL_CAPS inside constants.m is avoided; instead constants() returns a struct const.c, const.epsilon0, etc. [CODATA 2018].
  • Use MATLAB classdef, explicit property blocks, and argument validation (arguments blocks).

• Performance

  • Prefer vectorized operations.
  • Avoid per-pixel for loops in hot paths (FFT-based propagation scales as O(N log N) already).
  • For strongly diverging beams or large NA content, pad the field before propagation to reduce wrap-around.

• Numerical assumptions

  • The propagation engine currently assumes a spatially invariant refractive index n.
  • No non-paraxial vector diffraction, no anisotropic media, no graded index, no nonlinear effects.
  • Evanescent components are supported in angular spectrum mode and decay exponentially with distance.

---

License

License to be defined (research / academic use expected). Until a license is added, all rights are reserved.

---

Author and Contact

Optilab Project Lead Developer: nouni2 GitHub: https://github.com/Nouni2

---

Reference

[CODATA 2018] Mohr, P.J., Newell, D.B., Taylor, B.N. (2018). CODATA recommended values of the fundamental physical constants: 2018.