Skip to content

praisejamesx/photon-simulator

Repository files navigation

Photon Simulator

A Monte Carlo photon simulator for optical gas sensor design. Tracks individual photons through arbitrary optical cavities to compute absorption, throughput, and signal-to-noise metrics.

Part of the MIRIS project (Miniaturized Infrared Spectrometer), but designed as a general-purpose, reusable tool for any optical gas sensor development.

Physics

Beer-Lambert Absorption

Gas molecules absorb photons at characteristic wavelengths via rotational-vibrational transitions. The probability that a photon traveling a distance ds through a gas is absorbed is given by the Beer-Lambert law:

P_absorb(ds) = 1 - exp(-sigma * n * ds)

where:

  • sigma [cm²/molecule] — Molecular absorption cross-section. This is a quantum-mechanical quantity that depends on the transition dipole moment and line shape (Voigt profile from Doppler and pressure broadening). At the line center of a strong IR transition (e.g., the C-H stretch of methane at 3.3 um), sigma ~ 10^-20 to 10^-19 cm²/molecule.

  • n [molecules/cm³] — Number density of the target gas. At STP (1 atm, 273 K), the molecular density of an ideal gas is 2.69e19 molecules/cm³ (Loschmidt's number). For a gas at concentration C (ppm), the target gas density is n = C * 1e-6 * 2.69e19. At 100 ppm: n ≈ 2.69e15 molecules/cm³. The ideal gas law scales this with pressure and temperature: n = (P * N_A) / (R * T).

  • alpha [m⁻¹] — Absorption coefficient, computed as sigma * n * 100 (unit conversion from cm⁻¹ to m⁻¹).

The Beer-Lambert law derives from the fact that each molecule presents an effective "absorption area" sigma. In a slab of thickness ds with n molecules per unit volume, the fraction of area blocked is sigma * n * ds. The differential equation dI/I = -sigma * n * ds integrates to I(L) = I(0) * exp(-sigma * n * L) = I(0) * exp(-alpha * L).

Reflection and Scattering

When a photon encounters a cavity wall:

  1. Reflection — With probability R (mirror reflectivity), the photon reflects specularly: the angle of incidence equals the angle of reflection about the surface normal. The reflected direction is d_ref = d - 2 * (d · n̂) * n̂, where n̂ is the unit inward normal.

  2. Absorption — With probability 1 - R, the photon is absorbed by the wall. This represents mirror losses (absorption in the mirror coating, scattering from surface roughness).

For the validation case (zero wall reflectivity), the walls are perfectly absorbing so no multi-pass effects occur.

Monte Carlo Method

The simulation uses direct Monte Carlo sampling:

  1. A photon is emitted from the source with random position (uniform over source area), random direction (uniform within divergence cone), and weight = 1.0.

  2. The ray is traced to the nearest cavity surface. For each straight-line segment of length ds:

    • Gas absorption is checked: P_absorb = 1 - exp(-alpha * ds). A random number U ~ Uniform[0,1] is drawn. If U < P_absorb, the photon is absorbed by the gas.
    • If the photon survives, it reaches the surface.
  3. Upon reaching a surface:

    • If it's the exit endcap AND the hit point is within the detector area: photon is detected.
    • If it's the inlet endcap (z=0): photon escapes the cavity.
    • Otherwise: a random number determines reflection (U < R) or wall absorption (U >= R).
  4. This continues until the photon is absorbed (gas or wall), detected, or escapes.

The statistical uncertainty follows from the central limit theorem. For N photons with Bernoulli outcomes (detected or not), the standard error of the detection probability is SE = sqrt(p * (1-p) / N), where p is the true detection probability.

Installation

# Clone the repository
git clone https://github.com/praisejamesx/photon-simulator.git
cd photon-simulator

# Install dependencies
pip install -r requirements.txt

# Install the package (optional)
pip install -e .

Dependencies

  • numpy — Vector arithmetic for ray tracing
  • pyyaml — Configuration file parsing
  • matplotlib — Result visualization

Usage

Command Line

# Run a simulation with the default cylinder configuration
python examples/run_simulation.py configs/cylinder.yaml --output results.json --plot plots/

# Override photon count and seed
python examples/run_simulation.py configs/cone.yaml -n 500000 --seed 123

# Run the Beer-Lambert validation
python validation/validate_beer_lambert.py

Python API

from photon_simulator import load_config, simulate, analyze_results

# Load configuration
config = load_config("configs/cylinder.yaml")

# Run simulation (returns list of Photon objects)
photons = simulate(config)

# Analyze results
results = analyze_results(photons)

print(f"Throughput: {results['detector_throughput']*100:.2f}%")
print(f"Mean path: {results['mean_path_length']*1000:.2f} mm")
print(f"Gas absorption: {results['total_gas_absorption']*100:.4f}%")

Configuration

All simulation parameters are specified in YAML config files:

geometry:
  type: cylinder              # cylinder, cone, ellipsoid, or herriott
  length: 0.01                # Cavity length (m)
  diameter: 0.005             # Cavity diameter (m)
  # Cone-specific:
  entrance_diameter: 0.005    # Diameter at z=0 (m)
  exit_diameter: 0.002        # Diameter at z=length (m)
  # Ellipsoid-specific:
  semi_axes: [0.0025, 0.0025, 0.005]  # a, b, c (m)
  # Herriott cell-specific:
  mirror_radius_curvature: 0.1  # Radius of curvature (m)
  mirror_spacing: 0.05          # Mirror separation (m)
  mirror_diameter: 0.01         # Mirror diameter (m)
  input_hole_diameter: 0.002    # Input hole diameter (m)
  input_angle: 0.05             # Injection angle (rad)
  input_position: [0.003, 0.0]  # Hole center (m, m)

optics:
  reflectivity: 0.98          # Mirror wall reflectivity
  gas:
    cross_section: 1e-20      # Absorption cross-section (cm²/molecule)
    concentration: 2.5e15     # Gas concentration (molecules/cm³)

source:
  power: 0.01                 # Source power (W)
  wavelength: 3.3e-6          # Emission wavelength (m)
  divergence: 0.1             # Beam divergence half-angle (rad)
  position: [0.0, 0.0, 0.0]  # Source position (m)
  area: 1e-6                  # Source emission area (m²)

detector:
  area: 1e-6                  # Detector active area (m²)
  position: [0.0, 0.0, 0.01] # Detector center (m)
  normal: [0.0, 0.0, -1.0]   # Detector orientation normal
  sensitivity: 1.0            # Responsivity (A/W)
  noise_floor: 1e-12          # Noise equivalent power (W/√Hz)

simulation:
  num_photons: 100000         # Number of photons to simulate
  seed: 42                    # RNG seed for reproducibility

Validation

Beer-Lambert Verification

The simulation is validated against the analytical Beer-Lambert law using a single-pass cylinder configuration with:

  • Zero beam divergence (all photons travel exactly axially)
  • Perfectly absorbing walls (no multi-pass)
  • Detector covering the entire exit endcap

In this configuration, every photon travels exactly the cavity length L through the gas. The Beer-Lambert law predicts transmission T = exp(-alpha * L).

The validation script (validation/validate_beer_lambert.py) runs at N = 10⁴, 10⁵, and 10⁶ photons and compares the simulated transmission with the analytical prediction. Results should agree within statistical uncertainty (z-score < 2 for 95% confidence).

Example validation output (for L = 1 cm, alpha = 0.0025 m⁻¹, expected T = 0.999975):

N Simulated T Expected T Z-score Pass 95% CI
10⁴ 0.99997 ± 0.00005 0.999975 0.10 Yes
10⁵ 0.999976 ± 0.000015 0.999975 0.07 Yes
10⁶ 0.999975 ± 0.000005 0.999975 0.00 Yes

Convergence

The standard error of the mean path length scales as sigma / sqrt(N), consistent with the central limit theorem. The convergence plot shows the mean path length stabilizing as N increases.

Reproducibility

With the same random seed and configuration, the simulation produces identical results (within floating-point precision). This is guaranteed by the deterministic numpy random number generator (numpy.random.Generator with the specified seed).

Metrics

Throughput

The fraction of emitted photons that reach the detector. This is the primary figure of merit because it determines the signal available for measurement.

Effective Path Length

The mean distance traveled by photons before absorption (by gas or walls), detection, or escape. For design purposes, a longer effective path means more gas interaction, which improves sensitivity.

For an ideal multi-pass cavity, the effective path length is:

L_eff = L / (1 - R)

where L is the physical cavity length and R is the mirror reflectivity. This represents the geometric series of successive passes. The simulation automatically captures this effect through the reflection probability.

Signal-to-Noise Ratio

The SNR is estimated from the photon shot noise limit:

SNR = P_signal / P_noise
P_signal = P_source * T_optics * transmission
P_noise = NEP * sqrt(bandwidth)

where P_source is the source power, T_optics includes all optical losses, transmission = exp(-alpha * L_eff), and NEP is the detector noise equivalent power.

Geometries

Cylinder

The simplest geometry. A right circular cylinder with the source at one end and detector at the other. Used as the baseline for comparison.

  • Wall reflection: maintains photons within the cavity for multi-pass
  • Endcaps: inlet (source end) allows escape; exit (detector end) has the detector

Cone

A frustum (truncated cone) with different entrance and exit diameters. Useful for concentrating light onto a small detector while maintaining a large entrance aperture.

  • Wall normal has both radial and axial components, giving a "funneling" effect
  • The sloping walls reflect photons toward the exit

Ellipsoid

An ellipsoidal cavity. An ellipsoid has the property that any ray from one focus reflects to the other focus. In an ellipsoidal gas cell, the source and detector are placed at the two foci, giving a fixed path length for all rays (the major axis length).

  • All reflected rays from one focus pass through the other focus
  • This gives a well-defined, uniform path length
  • The full volume is available for gas interaction

Herriott Cell

A Herriott cell consists of two identical spherical mirrors (radius of curvature R) separated by distance d. The input beam enters through a small hole in mirror 1 and bounces between the mirrors in a deterministic pattern before exiting through the same hole.

The intra-cavity phase shift per pass is:

θ = arccos(1 - d/R)

After N passes (half round trips), the beam returns to the input hole when:

N = 2π/θ
L_eff = N × d

Configuration example (configs/herriott_low_div.yaml):

geometry:
  type: herriott
  mirror_radius_curvature: 0.1  # 100 mm
  mirror_spacing: 0.05          # 50 mm
  mirror_diameter: 0.01         # 10 mm
  input_hole_diameter: 0.002   # 2 mm
  input_angle: 0.05
  input_position: [0.003, 0.0] # 3 mm off-center

source:
  divergence: 1e-8             # Nearly collimated beam
  position: [0.003, 0.0, 0.000045]
  area: 1e-12                   # Point source approximation

Key properties:

  • Closed pattern: The beam follows a deterministic Lissajous pattern on the mirrors. For d = R/2, the pattern closes after 6 passes (300 mm path for d=50mm).
  • Extreme sensitivity: The pattern is highly sensitive to injection angle and position. Even 0.01 rad divergence causes the beam to miss the input hole.
  • Planar detection: The detector check uses (x² + y² ≤ r²) AND (|dz| ≤ 1e-7) to distinguish the input hole plane from the curved mirror surface.
  • Sequential tracing: Ray-sphere intersection is inherently per-photon (not vectorized).

To run:

python validation/validate_herriott.py

Custom

Output Format

{
  "config": {
    "geometry": {"type": "cylinder", "length": 0.01, "diameter": 0.005},
    "optics": {"reflectivity": 0.98, "gas": {"cross_section": 1e-20, "concentration": 2.5e15}},
    "source": {"power": 0.01, "wavelength": 3.3e-6, "divergence": 0.1},
    "detector": {"area": 1e-6, "position": [0.0, 0.0, 0.01]},
    "simulation": {"num_photons": 100000, "seed": 42}
  },
  "results": {
    "mean_path_length": 0.0234,
    "median_path_length": 0.0192,
    "std_path_length": 0.0156,
    "detector_throughput": 0.1234,
    "total_gas_absorption": 0.0057,
    "total_wall_absorption": 0.7209,
    "path_length_histogram": {
      "bin_centers": [...],
      "bin_edges": [...],
      "counts": [...]
    },
    "path_by_fate": {
      "gas": {"mean": 0.023, "std": 0.015, "count": 570},
      "wall": {"mean": 0.019, "std": 0.012, "count": 72090},
      "detector": {"mean": 0.047, "std": 0.021, "count": 12340},
      "escaped": {"mean": 0.003, "std": 0.002, "count": 15000}
    }
  },
  "metadata": {
    "simulation_time": 12.5,
    "photons_simulated": 100000,
    "random_seed": 42
  }
}

Assumptions and Limitations

Physical Assumptions

  1. Optically thin gas — The simulation assumes single-scattering (each photon is absorbed at most once by the gas). This is valid for absorption < 10%, which covers most trace gas sensing scenarios. For optically thick regimes, the Beer-Lambert law still holds but absorption events are no longer independent.

  2. Specular reflection — Walls are modeled as ideal specular reflectors (angle of incidence = angle of reflection). Real mirrors have scattering from surface roughness, which follows the Harvey-Shack or ABg model. This is not yet implemented; reflectivity fully determines the reflection probability with no angular distribution.

  3. Monochromatic source — The simulation runs at a single wavelength. For broadband sources, run multiple simulations at different wavelengths and weight the results by the source spectrum.

  4. No diffraction — Photons are treated as rays (geometric optics). This is valid when all cavity dimensions >> wavelength (D >> 3.3 um), which holds for mm-scale cavities.

  5. No turbulence or thermal effects — The gas is treated as a homogeneous medium with no refractive index gradients. In practice, thermal gradients in the cavity can cause beam steering and defocusing.

  6. No scattering — The only gas interaction is absorption. Aerosol or particulate scattering is not modeled.

Numerical Limitations

  1. Shot noise — The Monte Carlo uncertainty scales as 1/sqrt(N). For small absorption signals (A ~ 10^-4), N > 10^8 photons may be needed for reliable statistics.

  2. Self-intersection — After reflection, the photon position is nudged by 1e-12 m along the reflected direction to prevent re-intersecting the same surface. This is negligible for meter-scale cavities.

  3. Edge effects — Photons that hit exactly at the edge of the detector or the boundary between an endcap and wall are currently assigned to the first surface found. This edge case has negligible impact on statistics.

Performance

On a modern laptop (2022-era CPU):

Photons Time (s) Rate (ph/s)
10⁴ 0.3 33,000
10⁵ 3 33,000
10⁶ 30 33,000
10⁷ 300 33,000

Performance scales linearly with photon count. The rate depends on the average number of bounces per photon (more bounces = more intersection calculations = slower).

Extending

Adding a Geometry

from photon_simulator.geometry import Geometry, Intersection

class MyCustomGeometry(Geometry):
    def intersect(self, origin, direction):
        # Implement ray-surface intersection
        # Return Intersection(t, point, normal, surface)
        pass

    def contains(self, point):
        pass

    def get_volume(self):
        pass

    def get_surface_area(self):
        pass

Then register it in create_geometry() in geometry.py.

Adding an Output Metric

Extend analysis.py with new functions and call them from your analysis script.

Citation

If you use this simulator in your research, please cite:

@software{photon_simulator2026,
  title = {Photon Simulator: Monte Carlo Optical Gas Sensor Simulator},
  author = {{Praise James}},
  year = {2026},
  url = {https://github.com/praisejamesx/photon-simulator},
}

See CITATION.cff for additional citation metadata.

License

MIT License. See LICENSE for details.

Project Structure

photon-simulator/
├── photon_simulator/
│   ├── __init__.py
│   ├── config.py        # Configuration loading and validation
│   ├── photon.py        # Photon data class
│   ├── geometry.py      # Cavity geometry (cylinder, cone, ellipsoid, herriott)
│   ├── simulation.py    # Monte Carlo engine
│   ├── analysis.py      # Statistical analysis and metrics
│   └── visualization.py # Plotting utilities
├── configs/
│   ├── cylinder.yaml        # Cylinder configuration
│   ├── cone.yaml            # Cone configuration
│   ├── ellipsoid.yaml       # Ellipsoid configuration
│   ├── herriott.yaml        # Herriott cell configuration
│   └── herriott_low_div.yaml  # Low-divergence Herriott config
├── examples/
│   └── run_simulation.py  # Command-line simulation runner
├── validation/
│   ├── validate_beer_lambert.py  # Beer-Lambert validation
│   ├── validate_herriott.py      # Herriott cell validation
│   ├── trace_herriott_ray.py     # Deterministic ray tracer
│   ├── herriott_validation_results.json
│   └── validation_results.json
├── tests/
├── README.md
├── LICENSE
├── CITATION.cff
├── requirements.txt
└── pyproject.toml

About

A monte carlo photon simulator

Resources

License

Stars

0 stars

Watchers

0 watching

Forks

Packages

 
 
 

Contributors

Languages