# The Clear BSD License
#
# Copyright (c) 2026 Tobias Heibges
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted (subject to the limitations in the disclaimer
# below) provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from this
# software without specific prior written permission.
#
# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY
# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
"""
Gaussian Beam Source Implementation
Provides a Gaussian beam source following paraxial beam optics.
Suitable for modeling focused laser beams.
Examples
--------
>>> from surface_roughness.sources import GaussianBeam
>>>
>>> source = GaussianBeam(
... waist_position=(0, 0, 0),
... direction=(0, 0, 1),
... waist_radius=1e-3, # 1 mm waist
... num_rays=5000,
... wavelength=1064e-9, # Nd:YAG
... power=10e-3
... )
>>> rays = source.generate()
References
----------
.. [1] Siegman, A. E. (1986). Lasers. University Science Books.
"""
import numpy as np
from ..utilities.ray_data import RayBatch
from .base import RaySource
[docs]
class GaussianBeam(RaySource):
"""
Gaussian beam with specified waist and Rayleigh range.
Implements paraxial Gaussian beam using ray approximation.
Each ray represents a wavefront normal, with intensity weighted
according to the Gaussian profile.
Parameters
----------
waist_position : tuple of float
Position of beam waist (x, y, z) in meters.
direction : tuple of float
Beam axis direction (dx, dy, dz), will be normalized.
waist_radius : float
Beam waist radius (w₀) in meters.
num_rays : int
Number of rays to generate.
wavelength : float
Wavelength in meters. Must be monochromatic for Gaussian beam.
power : float, optional
Total beam power in watts. Default is 1.0.
Attributes
----------
waist_position : ndarray, shape (3,)
Beam waist position.
direction : ndarray, shape (3,)
Normalized beam direction.
waist_radius : float
Beam waist radius w₀.
rayleigh_range : float
Rayleigh range z_R = πw₀²/λ.
Notes
-----
The Gaussian beam has intensity profile:
I(r) = I₀ exp(-2r²/w²)
where w is the beam radius at distance z from the waist:
w(z) = w₀ √(1 + (z/z_R)²)
and z_R = πw₀²/λ is the Rayleigh range.
This implementation generates rays at the waist position with
parallel directions (plane wavefront at waist). Ray intensities
are weighted according to the Gaussian profile.
Examples
--------
>>> # 1 mm waist Nd:YAG laser
>>> source = GaussianBeam(
... waist_position=(0, 0, 0),
... direction=(0, 0, 1),
... waist_radius=1e-3,
... num_rays=5000,
... wavelength=1064e-9,
... power=10e-3
... )
>>> print(f"Rayleigh range: {source.rayleigh_range:.3f} m")
"""
[docs]
def __init__(
self,
waist_position: tuple[float, float, float],
direction: tuple[float, float, float],
waist_radius: float,
num_rays: int,
wavelength: float,
power: float = 1.0,
):
"""
Initialize Gaussian beam.
Parameters
----------
waist_position : tuple of float
Position of beam waist (x, y, z) in meters.
direction : tuple of float
Beam axis direction, will be normalized.
waist_radius : float
Beam waist radius w₀ in meters.
num_rays : int
Number of rays to generate.
wavelength : float
Wavelength in meters.
power : float, optional
Total beam power in watts. Default is 1.0.
Raises
------
ValueError
If wavelength is a tuple (polychromatic not supported),
or if waist_radius <= 0.
"""
if isinstance(wavelength, tuple):
raise ValueError("Gaussian beam requires monochromatic wavelength")
super().__init__(num_rays, wavelength, power)
self.waist_position = np.array(waist_position, dtype=np.float32)
self.waist_radius = waist_radius
# Normalize direction
direction_arr = np.array(direction, dtype=np.float32)
self.direction = direction_arr / np.linalg.norm(direction_arr)
# Compute Rayleigh range
self.rayleigh_range = np.pi * waist_radius**2 / wavelength
if waist_radius <= 0:
raise ValueError("waist_radius must be positive")
[docs]
def generate(self) -> RayBatch:
"""
Generate Gaussian beam.
Creates rays at the waist position with Gaussian-distributed
positions and parallel directions.
Returns
-------
RayBatch
Ray batch with Gaussian intensity distribution.
Notes
-----
Positions are sampled from a 2D Gaussian distribution with
σ = w₀/2. Intensities are weighted by the Gaussian profile
to accurately represent the beam's energy distribution.
"""
rays = self._allocate_rays()
# Create perpendicular basis
v1, v2 = self._create_perpendicular_basis(self.direction)
# Generate positions at waist with Gaussian distribution
# Use w0/2 as sigma for sampling to get good coverage
sigma = self.waist_radius / 2
x_local = np.random.normal(0, sigma, self.num_rays)
y_local = np.random.normal(0, sigma, self.num_rays)
# Positions at waist
rays.positions[:] = (
self.waist_position
+ x_local[:, np.newaxis] * v1
+ y_local[:, np.newaxis] * v2
)
# At waist, rays are parallel to axis (plane wavefront)
rays.directions[:] = self.direction
# Adjust intensities according to Gaussian profile
# I(r) ∝ exp(-2r²/w₀²)
r_squared = x_local**2 + y_local**2
rays.intensities[:] *= np.exp(-2 * r_squared / self.waist_radius**2)
# Renormalize to conserve power
rays.intensities[:] *= self.power / np.sum(rays.intensities)
# Monochromatic wavelength
rays.wavelengths[:] = self.wavelength
return rays
[docs]
def beam_radius_at(self, z: float) -> float:
"""
Compute beam radius at distance z from waist.
Parameters
----------
z : float
Distance from waist along beam axis in meters.
Returns
-------
float
Beam radius w(z) in meters.
Notes
-----
w(z) = w₀ √(1 + (z/z_R)²)
"""
return self.waist_radius * np.sqrt(1 + (z / self.rayleigh_range) ** 2)
[docs]
def divergence_angle(self) -> float:
"""
Compute far-field divergence angle.
Returns
-------
float
Half-angle divergence in radians.
Notes
-----
θ = λ / (π w₀)
"""
return self.wavelength / (np.pi * self.waist_radius)
[docs]
def __repr__(self) -> str:
"""Return string representation."""
return (
f"GaussianBeam(waist_position={self.waist_position.tolist()}, "
f"waist_radius={self.waist_radius:.2e}, "
f"rayleigh_range={self.rayleigh_range:.3f})"
)