psfsim.quadrature_integration
Gaussian quadrature integration optimized for exponential decay in detector.
Uses custom orthogonal polynomials built for the weight function w(z) = exp(-alpha*z) on the finite interval [0, detector_thickness], following the standard Gaussian quadrature construction method with three-term recurrence relations.
Classes
Build orthogonal polynomials optimized for exponential decay weight function. |
|
Adaptive Gaussian quadrature integrator for decaying functions. |
Functions
|
Build Gaussian quadrature nodes and weights for exponential decay weight function. |
Module Contents
- class ExponentialDecayPolynomials(alpha, z_min, z_max, n_order)[source]
Build orthogonal polynomials optimized for exponential decay weight function.
Uses Gram-Schmidt orthogonalization with the weight function w(z) = exp(-alpha*z) on the finite interval [z_min, z_max]. The resulting polynomials are used to construct Gaussian quadrature nodes and weights via the Golub-Welsch algorithm.
Computes the three-term recurrence relation coefficients and builds the Jacobi matrix for eigenvalue decomposition.
- Parameters:
alpha (float) – Decay constant in weight function w(z) = exp(-alpha*z).
z_min (float) – Integration interval bounds.
z_max (float) – Integration interval bounds.
n_order (int) – Maximum polynomial order to compute.
- _compute_moments()[source]
Compute moments m_k = integral of z^k * exp(-alpha*z) on [z_min, z_max].
Uses analytical formulas via integration by parts.
- Returns:
moments – moments[k] = integral of z^k * exp(-alpha*z) on [z_min, z_max]
- Return type:
np.ndarray of shape (2*n_order + 3,)
- _inner_product(p_coeffs, q_coeffs)[source]
Compute inner product of two polynomials with exponential weight.
<p, q> = integral of p(z) * q(z) * exp(-alpha*z) dz
Uses moment-based formula for efficiency: <p, q> = sum_i sum_j p_i * q_j * m_{i+j}
Uses einsum for vectorized computation.
- Parameters:
p_coeffs (np.ndarray) – Coefficients of polynomial p in ascending order of powers.
q_coeffs (np.ndarray) – Coefficients of polynomial q in ascending order of powers.
- Returns:
ip – Inner product.
- Return type:
float
- _monic_poly_multiply_z(p_asc)[source]
Multiply polynomial by z: if p = sum_i p_i z^i, return sum_i p_i z^{i+1}.
- Parameters:
p_asc (np.ndarray) – Coefficients in ascending order of powers.
- Returns:
result – Coefficients of z*p in ascending order.
- Return type:
np.ndarray
- _build_recurrence_coefficients()[source]
Build three-term recurrence relation coefficients.
Standard form: p_{n+1}(z) = (z - alpha_n) * p_n(z) - beta_n * p_{n-1}(z)
- where:
alpha_n = <z*p_n, p_n> / <p_n, p_n> beta_n = <p_n, p_n> / <p_{n-1}, p_{n-1}>
This builds the Jacobi matrix diagonals for the Golub-Welsch algorithm.
- compute_quadrature_nodes_and_weights(n)[source]
Compute Gaussian quadrature nodes and weights using Golub-Welsch algorithm.
The Golub-Welsch algorithm:
Build the symmetric tridiagonal Jacobi matrix from recurrence coefficients
Compute eigenvalues (which are the quadrature nodes)
Compute weights from the first row of the eigenvector matrix
- Parameters:
n (int) – Number of quadrature points.
- Returns:
nodes (np.ndarray of shape (n,)) – Quadrature nodes (sorted, in [z_min, z_max]).
weights (np.ndarray of shape (n,)) – Quadrature weights (corresponding to nodes). These are in the form of int_a^b f(z) dz = sum_j v_j f(z_j).
- build_exponential_decay_quadrature(alpha, z_min, z_max, n_order)[source]
Build Gaussian quadrature nodes and weights for exponential decay weight function.
Constructs orthogonal polynomials with respect to the weight function w(z) = exp(-alpha*z) on the finite interval [z_min, z_max], then extracts quadrature nodes and weights using the Golub-Welsch algorithm.
- Parameters:
alpha (float) – Decay constant in weight function w(z) = exp(-alpha*z).
z_min (float) – Integration interval.
z_max (float) – Integration interval.
n_order (int) – Number of quadrature points.
- Returns:
nodes (np.ndarray of shape (n_order,)) – Quadrature nodes in [z_min, z_max].
weights (np.ndarray of shape (n_order,)) – Quadrature weights.
- class QuadratureIntegrator(wavelength, detector_thickness, ux, uy, filter_obj)[source]
Adaptive Gaussian quadrature integrator for decaying functions.
Optimizes quadrature nodes and weights based on the characteristic decay length of the intensity function in the detector. Uses custom orthogonal polynomials built for exponential decay weight function.
- ux, uy
Orthographic projections of ray directions (normalized pupil coordinates).
- Type:
np.ndarray
- _compute_kz_imag()[source]
Compute the imaginary part of kz (decay constant) across spatial grid.
- Returns:
kz_imag – Imaginary part of wave vector in detector material. Positive values indicate exponential decay: I(z) ~ exp(-2*kz_imag*z).
- Return type:
np.ndarray (same shape as ux, uy)
- analyze_decay()[source]
Analyze decay behavior of intensity across spatial grid.
- Returns:
decay_length_characteristic (float) – Characteristic decay length (1/alpha) in microns.
decay_constant (float) – Characteristic decay constant alpha in 1/microns.
- _adaptive_order(decay_length)[source]
Select quadrature order based on decay length relative to detector thickness.
- Parameters:
decay_length (float) – Characteristic decay length in microns.
- Returns:
order – Number of quadrature points (3 to 20).
- Return type:
int
- get_nodes_and_weights()[source]
Get optimized quadrature nodes and weights.
Uses custom orthogonal polynomials built for the exponential decay weight function w(z) = exp(-alpha*z) on [0, detector_thickness].
- Returns:
z_nodes (np.ndarray) – Depth coordinates (in microns) for quadrature points.
weights (np.ndarray) – Quadrature weights.
order (int) – Number of quadrature points used.
- integrate(intensity_array, axis=2)[source]
Integrate intensity along specified axis using custom quadrature.
- Parameters:
intensity_array (np.ndarray) – Intensity values with shape (…, nz) where nz is the number of z-points.
axis (int, optional) – Axis along which to integrate (default: 2, for shape (ux, uy, z)).
- Returns:
integrated – Integrated intensity with shape matching intensity_array without the integration axis.
- Return type:
np.ndarray