# -*- encoding: utf-8 -*-
"""
Gaussian optics.
The module implements:
- Ray transfer matrices for geometrical and gaussian optics.
See RayTransferMatrix, GeometricRay and BeamParameter
- Conjugation relations for geometrical and gaussian optics.
See geometric_conj*, gauss_conj and conjugate_gauss_beams
The conventions for the distances are as follows:
focal distance
positive for convergent lenses
object distance
positive for real objects
image distance
positive for real images
"""
from __future__ import print_function, division
__all__ = [
'RayTransferMatrix',
'FreeSpace',
'FlatRefraction',
'CurvedRefraction',
'FlatMirror',
'CurvedMirror',
'ThinLens',
'GeometricRay',
'BeamParameter',
'waist2rayleigh',
'rayleigh2waist',
'geometric_conj_ab',
'geometric_conj_af',
'geometric_conj_bf',
'gaussian_conj',
'conjugate_gauss_beams',
]
from sympy import (atan2, Expr, I, im, Matrix, oo, pi, re, sqrt, sympify,
together)
from sympy.utilities.misc import filldedent
###
# A, B, C, D matrices
###
[docs]class RayTransferMatrix(Matrix):
"""
Base class for a Ray Transfer Matrix.
It should be used if there isn't already a more specific subclass mentioned
in See Also.
Parameters
==========
parameters : A, B, C and D or 2x2 matrix (Matrix(2, 2, [A, B, C, D]))
Examples
========
>>> from sympy.physics.optics import RayTransferMatrix, ThinLens
>>> from sympy import Symbol, Matrix
>>> mat = RayTransferMatrix(1, 2, 3, 4)
>>> mat
Matrix([
[1, 2],
[3, 4]])
>>> RayTransferMatrix(Matrix([[1, 2], [3, 4]]))
Matrix([
[1, 2],
[3, 4]])
>>> mat.A
1
>>> f = Symbol('f')
>>> lens = ThinLens(f)
>>> lens
Matrix([
[ 1, 0],
[-1/f, 1]])
>>> lens.C
-1/f
See Also
========
GeometricRay, BeamParameter,
FreeSpace, FlatRefraction, CurvedRefraction,
FlatMirror, CurvedMirror, ThinLens
References
==========
.. [1] http://en.wikipedia.org/wiki/Ray_transfer_matrix_analysis
"""
def __new__(cls, *args):
if len(args) == 4:
temp = ((args[0], args[1]), (args[2], args[3]))
elif len(args) == 1 \
and isinstance(args[0], Matrix) \
and args[0].shape == (2, 2):
temp = args[0]
else:
raise ValueError(filldedent('''
Expecting 2x2 Matrix or the 4 elements of
the Matrix but got %s''' % str(args)))
return Matrix.__new__(cls, temp)
def __mul__(self, other):
if isinstance(other, RayTransferMatrix):
return RayTransferMatrix(Matrix.__mul__(self, other))
elif isinstance(other, GeometricRay):
return GeometricRay(Matrix.__mul__(self, other))
elif isinstance(other, BeamParameter):
temp = self*Matrix(((other.q,), (1,)))
q = (temp[0]/temp[1]).expand(complex=True)
return BeamParameter(other.wavelen,
together(re(q)),
z_r=together(im(q)))
else:
return Matrix.__mul__(self, other)
@property
def A(self):
"""
The A parameter of the Matrix.
Examples
========
>>> from sympy.physics.optics import RayTransferMatrix
>>> mat = RayTransferMatrix(1, 2, 3, 4)
>>> mat.A
1
"""
return self[0, 0]
@property
def B(self):
"""
The B parameter of the Matrix.
Examples
========
>>> from sympy.physics.optics import RayTransferMatrix
>>> mat = RayTransferMatrix(1, 2, 3, 4)
>>> mat.B
2
"""
return self[0, 1]
@property
def C(self):
"""
The C parameter of the Matrix.
Examples
========
>>> from sympy.physics.optics import RayTransferMatrix
>>> mat = RayTransferMatrix(1, 2, 3, 4)
>>> mat.C
3
"""
return self[1, 0]
@property
def D(self):
"""
The D parameter of the Matrix.
Examples
========
>>> from sympy.physics.optics import RayTransferMatrix
>>> mat = RayTransferMatrix(1, 2, 3, 4)
>>> mat.D
4
"""
return self[1, 1]
[docs]class FreeSpace(RayTransferMatrix):
"""
Ray Transfer Matrix for free space.
Parameters
==========
distance
See Also
========
RayTransferMatrix
Examples
========
>>> from sympy.physics.optics import FreeSpace
>>> from sympy import symbols
>>> d = symbols('d')
>>> FreeSpace(d)
Matrix([
[1, d],
[0, 1]])
"""
def __new__(cls, d):
return RayTransferMatrix.__new__(cls, 1, d, 0, 1)
[docs]class FlatRefraction(RayTransferMatrix):
"""
Ray Transfer Matrix for refraction.
Parameters
==========
n1 : refractive index of one medium
n2 : refractive index of other medium
See Also
========
RayTransferMatrix
Examples
========
>>> from sympy.physics.optics import FlatRefraction
>>> from sympy import symbols
>>> n1, n2 = symbols('n1 n2')
>>> FlatRefraction(n1, n2)
Matrix([
[1, 0],
[0, n1/n2]])
"""
def __new__(cls, n1, n2):
n1, n2 = map(sympify, (n1, n2))
return RayTransferMatrix.__new__(cls, 1, 0, 0, n1/n2)
[docs]class CurvedRefraction(RayTransferMatrix):
"""
Ray Transfer Matrix for refraction on curved interface.
Parameters
==========
R : radius of curvature (positive for concave)
n1 : refractive index of one medium
n2 : refractive index of other medium
See Also
========
RayTransferMatrix
Examples
========
>>> from sympy.physics.optics import CurvedRefraction
>>> from sympy import symbols
>>> R, n1, n2 = symbols('R n1 n2')
>>> CurvedRefraction(R, n1, n2)
Matrix([
[ 1, 0],
[(n1 - n2)/(R*n2), n1/n2]])
"""
def __new__(cls, R, n1, n2):
R, n1, n2 = map(sympify, (R, n1, n2))
return RayTransferMatrix.__new__(cls, 1, 0, (n1 - n2)/R/n2, n1/n2)
[docs]class FlatMirror(RayTransferMatrix):
"""
Ray Transfer Matrix for reflection.
See Also
========
RayTransferMatrix
Examples
========
>>> from sympy.physics.optics import FlatMirror
>>> FlatMirror()
Matrix([
[1, 0],
[0, 1]])
"""
def __new__(cls):
return RayTransferMatrix.__new__(cls, 1, 0, 0, 1)
[docs]class CurvedMirror(RayTransferMatrix):
"""
Ray Transfer Matrix for reflection from curved surface.
Parameters
==========
R : radius of curvature (positive for concave)
See Also
========
RayTransferMatrix
Examples
========
>>> from sympy.physics.optics import CurvedMirror
>>> from sympy import symbols
>>> R = symbols('R')
>>> CurvedMirror(R)
Matrix([
[ 1, 0],
[-2/R, 1]])
"""
def __new__(cls, R):
R = sympify(R)
return RayTransferMatrix.__new__(cls, 1, 0, -2/R, 1)
[docs]class ThinLens(RayTransferMatrix):
"""
Ray Transfer Matrix for a thin lens.
Parameters
==========
f : the focal distance
See Also
========
RayTransferMatrix
Examples
========
>>> from sympy.physics.optics import ThinLens
>>> from sympy import symbols
>>> f = symbols('f')
>>> ThinLens(f)
Matrix([
[ 1, 0],
[-1/f, 1]])
"""
def __new__(cls, f):
f = sympify(f)
return RayTransferMatrix.__new__(cls, 1, 0, -1/f, 1)
###
# Representation for geometric ray
###
[docs]class GeometricRay(Matrix):
"""
Representation for a geometric ray in the Ray Transfer Matrix formalism.
Parameters
==========
h : height, and
angle : angle, or
matrix : a 2x1 matrix (Matrix(2, 1, [height, angle]))
Examples
========
>>> from sympy.physics.optics import GeometricRay, FreeSpace
>>> from sympy import symbols, Matrix
>>> d, h, angle = symbols('d, h, angle')
>>> GeometricRay(h, angle)
Matrix([
[ h],
[angle]])
>>> FreeSpace(d)*GeometricRay(h, angle)
Matrix([
[angle*d + h],
[ angle]])
>>> GeometricRay( Matrix( ((h,), (angle,)) ) )
Matrix([
[ h],
[angle]])
See Also
========
RayTransferMatrix
"""
def __new__(cls, *args):
if len(args) == 1 and isinstance(args[0], Matrix) \
and args[0].shape == (2, 1):
temp = args[0]
elif len(args) == 2:
temp = ((args[0],), (args[1],))
else:
raise ValueError(filldedent('''
Expecting 2x1 Matrix or the 2 elements of
the Matrix but got %s''' % str(args)))
return Matrix.__new__(cls, temp)
@property
def height(self):
"""
The distance from the optical axis.
Examples
========
>>> from sympy.physics.optics import GeometricRay
>>> from sympy import symbols
>>> h, angle = symbols('h, angle')
>>> gRay = GeometricRay(h, angle)
>>> gRay.height
h
"""
return self[0]
@property
def angle(self):
"""
The angle with the optical axis.
Examples
========
>>> from sympy.physics.optics import GeometricRay
>>> from sympy import symbols
>>> h, angle = symbols('h, angle')
>>> gRay = GeometricRay(h, angle)
>>> gRay.angle
angle
"""
return self[1]
###
# Representation for gauss beam
###
[docs]class BeamParameter(Expr):
"""
Representation for a gaussian ray in the Ray Transfer Matrix formalism.
Parameters
==========
wavelen : the wavelength,
z : the distance to waist, and
w : the waist, or
z_r : the rayleigh range
Examples
========
>>> from sympy.physics.optics import BeamParameter
>>> p = BeamParameter(530e-9, 1, w=1e-3)
>>> p.q
1 + 1.88679245283019*I*pi
>>> p.q.n()
1.0 + 5.92753330865999*I
>>> p.w_0.n()
0.00100000000000000
>>> p.z_r.n()
5.92753330865999
>>> from sympy.physics.optics import FreeSpace
>>> fs = FreeSpace(10)
>>> p1 = fs*p
>>> p.w.n()
0.00101413072159615
>>> p1.w.n()
0.00210803120913829
See Also
========
RayTransferMatrix
References
==========
.. [1] http://en.wikipedia.org/wiki/Complex_beam_parameter
.. [2] http://en.wikipedia.org/wiki/Gaussian_beam
"""
#TODO A class Complex may be implemented. The BeamParameter may
# subclass it. See:
# https://groups.google.com/d/topic/sympy/7XkU07NRBEs/discussion
__slots__ = ['z', 'z_r', 'wavelen']
def __new__(cls, wavelen, z, **kwargs):
wavelen, z = map(sympify, (wavelen, z))
inst = Expr.__new__(cls, wavelen, z)
inst.wavelen = wavelen
inst.z = z
if len(kwargs) != 1:
raise ValueError('Constructor expects exactly one named argument.')
elif 'z_r' in kwargs:
inst.z_r = sympify(kwargs['z_r'])
elif 'w' in kwargs:
inst.z_r = waist2rayleigh(sympify(kwargs['w']), wavelen)
else:
raise ValueError('The constructor needs named argument w or z_r')
return inst
@property
def q(self):
"""
The complex parameter representing the beam.
Examples
========
>>> from sympy.physics.optics import BeamParameter
>>> p = BeamParameter(530e-9, 1, w=1e-3)
>>> p.q
1 + 1.88679245283019*I*pi
"""
return self.z + I*self.z_r
@property
def radius(self):
"""
The radius of curvature of the phase front.
Examples
========
>>> from sympy.physics.optics import BeamParameter
>>> p = BeamParameter(530e-9, 1, w=1e-3)
>>> p.radius
1 + 3.55998576005696*pi**2
"""
return self.z*(1 + (self.z_r/self.z)**2)
@property
def w(self):
"""
The beam radius at `1/e^2` intensity.
See Also
========
w_0 : the minimal radius of beam
Examples
========
>>> from sympy.physics.optics import BeamParameter
>>> p = BeamParameter(530e-9, 1, w=1e-3)
>>> p.w
0.001*sqrt(0.2809/pi**2 + 1)
"""
return self.w_0*sqrt(1 + (self.z/self.z_r)**2)
@property
def w_0(self):
"""
The beam waist (minimal radius).
See Also
========
w : the beam radius at `1/e^2` intensity
Examples
========
>>> from sympy.physics.optics import BeamParameter
>>> p = BeamParameter(530e-9, 1, w=1e-3)
>>> p.w_0
0.00100000000000000
"""
return sqrt(self.z_r/pi*self.wavelen)
@property
def divergence(self):
"""
Half of the total angular spread.
Examples
========
>>> from sympy.physics.optics import BeamParameter
>>> p = BeamParameter(530e-9, 1, w=1e-3)
>>> p.divergence
0.00053/pi
"""
return self.wavelen/pi/self.w_0
@property
def gouy(self):
"""
The Gouy phase.
Examples
========
>>> from sympy.physics.optics import BeamParameter
>>> p = BeamParameter(530e-9, 1, w=1e-3)
>>> p.gouy
atan(0.53/pi)
"""
return atan2(self.z, self.z_r)
@property
def waist_approximation_limit(self):
"""
The minimal waist for which the gauss beam approximation is valid.
The gauss beam is a solution to the paraxial equation. For curvatures
that are too great it is not a valid approximation.
Examples
========
>>> from sympy.physics.optics import BeamParameter
>>> p = BeamParameter(530e-9, 1, w=1e-3)
>>> p.waist_approximation_limit
1.06e-6/pi
"""
return 2*self.wavelen/pi
###
# Utilities
###
[docs]def waist2rayleigh(w, wavelen):
"""
Calculate the rayleigh range from the waist of a gaussian beam.
See Also
========
rayleigh2waist, BeamParameter
Examples
========
>>> from sympy.physics.optics import waist2rayleigh
>>> from sympy import symbols
>>> w, wavelen = symbols('w wavelen')
>>> waist2rayleigh(w, wavelen)
pi*w**2/wavelen
"""
w, wavelen = map(sympify, (w, wavelen))
return w**2*pi/wavelen
[docs]def rayleigh2waist(z_r, wavelen):
"""Calculate the waist from the rayleigh range of a gaussian beam.
See Also
========
waist2rayleigh, BeamParameter
Examples
========
>>> from sympy.physics.optics import rayleigh2waist
>>> from sympy import symbols
>>> z_r, wavelen = symbols('z_r wavelen')
>>> rayleigh2waist(z_r, wavelen)
sqrt(wavelen*z_r)/sqrt(pi)
"""
z_r, wavelen = map(sympify, (z_r, wavelen))
return sqrt(z_r/pi*wavelen)
[docs]def geometric_conj_ab(a, b):
"""
Conjugation relation for geometrical beams under paraxial conditions.
Takes the distances to the optical element and returns the needed
focal distance.
See Also
========
geometric_conj_af, geometric_conj_bf
Examples
========
>>> from sympy.physics.optics import geometric_conj_ab
>>> from sympy import symbols
>>> a, b = symbols('a b')
>>> geometric_conj_ab(a, b)
a*b/(a + b)
"""
a, b = map(sympify, (a, b))
if abs(a) == oo or abs(b) == oo:
return a if abs(b) == oo else b
else:
return a*b/(a + b)
[docs]def geometric_conj_af(a, f):
"""
Conjugation relation for geometrical beams under paraxial conditions.
Takes the object distance (for geometric_conj_af) or the image distance
(for geometric_conj_bf) to the optical element and the focal distance.
Then it returns the other distance needed for conjugation.
See Also
========
geometric_conj_ab
Examples
========
>>> from sympy.physics.optics.gaussopt import geometric_conj_af, geometric_conj_bf
>>> from sympy import symbols
>>> a, b, f = symbols('a b f')
>>> geometric_conj_af(a, f)
a*f/(a - f)
>>> geometric_conj_bf(b, f)
b*f/(b - f)
"""
a, f = map(sympify, (a, f))
return -geometric_conj_ab(a, -f)
geometric_conj_bf = geometric_conj_af
[docs]def gaussian_conj(s_in, z_r_in, f):
"""
Conjugation relation for gaussian beams.
Parameters
==========
s_in : the distance to optical element from the waist
z_r_in : the rayleigh range of the incident beam
f : the focal length of the optical element
Returns
=======
a tuple containing (s_out, z_r_out, m)
s_out : the distance between the new waist and the optical element
z_r_out : the rayleigh range of the emergent beam
m : the ration between the new and the old waists
Examples
========
>>> from sympy.physics.optics import gaussian_conj
>>> from sympy import symbols
>>> s_in, z_r_in, f = symbols('s_in z_r_in f')
>>> gaussian_conj(s_in, z_r_in, f)[0]
1/(-1/(s_in + z_r_in**2/(-f + s_in)) + 1/f)
>>> gaussian_conj(s_in, z_r_in, f)[1]
z_r_in/(1 - s_in**2/f**2 + z_r_in**2/f**2)
>>> gaussian_conj(s_in, z_r_in, f)[2]
1/sqrt(1 - s_in**2/f**2 + z_r_in**2/f**2)
"""
s_in, z_r_in, f = map(sympify, (s_in, z_r_in, f))
s_out = 1 / ( -1/(s_in + z_r_in**2/(s_in - f)) + 1/f )
m = 1/sqrt((1 - (s_in/f)**2) + (z_r_in/f)**2)
z_r_out = z_r_in / ((1 - (s_in/f)**2) + (z_r_in/f)**2)
return (s_out, z_r_out, m)
[docs]def conjugate_gauss_beams(wavelen, waist_in, waist_out, **kwargs):
"""
Find the optical setup conjugating the object/image waists.
Parameters
==========
wavelen : the wavelength of the beam
waist_in and waist_out : the waists to be conjugated
f : the focal distance of the element used in the conjugation
Returns
=======
a tuple containing (s_in, s_out, f)
s_in : the distance before the optical element
s_out : the distance after the optical element
f : the focal distance of the optical element
Examples
========
>>> from sympy.physics.optics import conjugate_gauss_beams
>>> from sympy import symbols, factor
>>> l, w_i, w_o, f = symbols('l w_i w_o f')
>>> conjugate_gauss_beams(l, w_i, w_o, f=f)[0]
f*(-sqrt(w_i**2/w_o**2 - pi**2*w_i**4/(f**2*l**2)) + 1)
>>> factor(conjugate_gauss_beams(l, w_i, w_o, f=f)[1])
f*w_o**2*(w_i**2/w_o**2 - sqrt(w_i**2/w_o**2 -
pi**2*w_i**4/(f**2*l**2)))/w_i**2
>>> conjugate_gauss_beams(l, w_i, w_o, f=f)[2]
f
"""
#TODO add the other possible arguments
wavelen, waist_in, waist_out = map(sympify, (wavelen, waist_in, waist_out))
m = waist_out / waist_in
z = waist2rayleigh(waist_in, wavelen)
if len(kwargs) != 1:
raise ValueError("The function expects only one named argument")
elif 'dist' in kwargs:
raise NotImplementedError(filldedent('''
Currently only focal length is supported as a parameter'''))
elif 'f' in kwargs:
f = sympify(kwargs['f'])
s_in = f * (1 - sqrt(1/m**2 - z**2/f**2))
s_out = gaussian_conj(s_in, z, f)[0]
elif 's_in' in kwargs:
raise NotImplementedError(filldedent('''
Currently only focal length is supported as a parameter'''))
else:
raise ValueError(filldedent('''
The functions expects the focal length as a named argument'''))
return (s_in, s_out, f)
#TODO
#def plot_beam():
# """Plot the beam radius as it propagates in space."""
# pass
#TODO
#def plot_beam_conjugation():
# """
# Plot the intersection of two beams.
#
# Represents the conjugation relation.
#
# See Also
# ========
#
# conjugate_gauss_beams
# """
# pass
