deep_causality_physics/photonics/
ray.rs

1/*
2 * SPDX-License-Identifier: MIT
3 * Copyright (c) 2023 - 2026. The DeepCausality Authors and Contributors. All Rights Reserved.
4 */
5
6use crate::photonics::quantities::{AbcdMatrix, OpticalPower, RayAngle, RayHeight};
7
8use crate::{IndexOfRefraction, PhysicsError};
9use deep_causality_tensor::{CausalTensor, Tensor};
10
11/// Applies an ABCD matrix to a ray vector.
12///
13/// $$ \begin{pmatrix} y_{out} \\ \theta_{out} \end{pmatrix} = \begin{pmatrix} A & B \\ C & D \end{pmatrix} \begin{pmatrix} y_{in} \\ \theta_{in} \end{pmatrix} $$
14///
15/// The ray transfer matrix (ABCD matrix) describes the optical system in the paraxial approximation.
16/// The vector consists of the ray height $y$ and angle $\theta$.
17///
18/// # Arguments
19/// *   `matrix` - The $2 \times 2$ ray transfer matrix (ABCD).
20/// *   `height` - Input ray height $y$.
21/// *   `angle` - Input ray angle $\theta$.
22///
23/// # Returns
24/// *   `Result<(RayHeight, RayAngle), PhysicsError>` - The output ray height and angle.
25pub fn ray_transfer_kernel(
26    matrix: &AbcdMatrix,
27    height: RayHeight,
28    angle: RayAngle,
29) -> Result<(RayHeight, RayAngle), PhysicsError> {
30    let m = matrix.inner();
31    if m.shape() != [2, 2] {
32        return Err(PhysicsError::DimensionMismatch(
33            "ABCD matrix must be 2x2".into(),
34        ));
35    }
36
37    let y_in = height.value();
38    let theta_in = angle.value();
39
40    // Create input column vector [y, theta]
41    let input_vec = CausalTensor::new(vec![y_in, theta_in], vec![2, 1])?;
42
43    // Matrix multiplication: [y_out, theta_out] = M * [y_in, theta_in]
44    let output_vec = m.matmul(&input_vec)?; // output_vec will be [2, 1]
45
46    let y_out = output_vec.data()[0];
47    let theta_out = output_vec.data()[1];
48
49    Ok((RayHeight::new(y_out)?, RayAngle::new(theta_out)?))
50}
51
52/// Calculates refracted angle using Snell's Law or returns a Critical Angle error.
53///
54/// $$ n_1 \sin \theta_1 = n_2 \sin \theta_2 $$
55///
56/// # Arguments
57/// *   `n1` - Refractive index of medium 1.
58/// *   `n2` - Refractive index of medium 2.
59/// *   `theta1` - Angle of incidence (relative to normal).
60///
61/// # Returns
62/// *   `Result<RayAngle, PhysicsError>` - Angle of refraction.
63///
64/// # Errors
65/// *   `PhysicalInvariantBroken` - If total internal reflection occurs ($\sin \theta_2 > 1$).
66pub fn snells_law_kernel(
67    n1: IndexOfRefraction,
68    n2: IndexOfRefraction,
69    theta1: RayAngle,
70) -> Result<RayAngle, PhysicsError> {
71    let n1_val = n1.value();
72    let n2_val = n2.value();
73    let theta1_val = theta1.value();
74
75    // sin(theta2) = (n1 / n2) * sin(theta1)
76    let sin_theta2 = (n1_val / n2_val) * theta1_val.sin();
77
78    if sin_theta2.abs() > 1.0 {
79        return Err(PhysicsError::PhysicalInvariantBroken(
80            "Total Internal Reflection: sin(theta2) > 1".into(),
81        ));
82    }
83
84    let theta2 = sin_theta2.asin();
85    RayAngle::new(theta2)
86}
87
88/// Calculates optical power and focal length using the Lens Maker's Equation.
89///
90/// $$ P = (n - 1) \left( \frac{1}{R_1} - \frac{1}{R_2} \right) $$
91///
92/// Uses the sign convention where:
93/// *   Light travels from left to right.
94/// *   Radius of curvature $R$ is positive if the center of curvature is to the right of the surface.
95/// *   Therefore, for a biconvex lens, $R_1 > 0$ (front surface convex) and $R_2 < 0$ (back surface convex).
96///
97/// # Arguments
98/// *   `n` - Refractive index of the lens material (assuming ambient index is 1).
99/// *   `r1` - Radius of curvature of the first surface.
100/// *   `r2` - Radius of curvature of the second surface.
101///
102/// # Returns
103/// *   `Result<OpticalPower, PhysicsError>` - Optical power in Diopters.
104pub fn lens_maker_kernel(
105    n: IndexOfRefraction,
106    r1: f64,
107    r2: f64,
108) -> Result<OpticalPower, PhysicsError> {
109    // P = (n - 1) * (1/R1 - 1/R2)
110    let n_val = n.value();
111
112    if r1 == 0.0 || r2 == 0.0 {
113        return Err(PhysicsError::Singularity(
114            "Radius of curvature cannot be zero".into(),
115        ));
116    }
117
118    let power = (n_val - 1.0) * ((1.0 / r1) - (1.0 / r2));
119    OpticalPower::new(power)
120}
deep_causality_physics/photonics/ray.rs

deep_causality_physics/photonics/
ray.rs
