Repository: DCC-Lab/RayTracing

Issue:
Axicon does not change the values of backIndex and frontIndex
The `Axicon` class does not change the default values of `backIndex` and `frontIndex` with the input value `n`.

Discussion:

Patch:
diff --git a/raytracing/axicon.py b/raytracing/axicon.py
index 4bf6013..680a686 100644
--- a/raytracing/axicon.py
+++ b/raytracing/axicon.py
@@ -2,130 +2,170 @@ from .matrix import *
 import matplotlib.pyplot as plt
 
 
-
 class Axicon(Matrix):
-	"""
-	This class is an advanced module that describes an axicon lens, not part of the basic formalism.
-	Using this class an axicon conical lens can be presented.
-	Axicon lenses are used to obtain a line focus instead of a point.
-
-	Parameters
-	----------
-	alpha : float
-		alpha is the small angle in radians of the axicon
-		(typically 2.5 or 5 degrees) corresponding to 90-apex angle
-	n : float
-		index of refraction.
-		This value cannot be less than 1.0.
-	diameter : float
-		Aperture of the element. (default = +Inf)
-		The diameter of the aperture must be a positive value.
-	label : string
-		The label of the axicon lens.
-
-	"""
-
-	def __init__(self, alpha, n, diameter=float('+Inf'), label=''):	
-
-		self.n = n
-		self.alpha = alpha
-		super(Axicon, self).__init__(A=1, B=0, C=0,D=1, physicalLength=0, apertureDiameter=diameter, label=label)
-
-	def deviationAngle(self):
-		""" This function provides deviation angle delta
-
-		Returns
-		-------
-		delta : float
-			the deviation angle
-
-		See ALso
-		--------
-		https://ru.b-ok2.org/book/2482970/9062b7, p.48
-
-		"""
-
-		return (self.n-1.0)*self.alpha
-
-
-	def focalLineLength(self, yMax=None):
-		""" Provides the line length, assuming a ray at height yMax
-
-		Parameters
-		----------
-		yMax : float
-			the height of the ray (default=None)
-			If no height is defined for the ray, then yMax would be set to the height of the axicon (apertureDiameter/2)
-
-		Returns
-		-------
-		focalLineLength : float
-			the length of the focal line
-
-		See ALso
-		--------
-		https://ru.b-ok2.org/book/2482970/9062b7, p.48
-
-		"""
+    """
+    This class is an advanced module that describes an axicon lens, 
+    not part of the basic formalism. Using this class an axicon 
+    conical lens can be presented.
+    Axicon lenses are used to obtain a line focus instead of a point.
+    The formalism is described in Kloos, sec. 2.2.3
 
-		if yMax == None:
-			yMax = self.apertureDiameter/2
+    Parameters
+    ----------
+    alpha : float
+        alpha is the small angle in radians of the axicon
+        (typically 2.5 or 5 degrees) corresponding to 90-apex angle
+    n : float
+        index of refraction.
+        This value cannot be less than 1.0. It is assumed the axicon
+        is in air.
+    diameter : float
+        Aperture of the element. (default = +Inf)
+        The diameter of the aperture must be a positive value.
+    label : string
+        The label of the axicon lens.
 
-		return yMax/(self.n-1.0)/self.alpha
+    """
 
-	def mul_ray(self, rightSideRay):
-		""" This function is used to calculate the output ray through an axicon.
+    def __init__(self, alpha, n, diameter=float('+Inf'), label=''):
 
-		Parameters
-		----------
-		rightSideRay : object of ray class
-			A ray with a defined height and angle.
+        self.n = n
+        self.alpha = alpha
+        super(Axicon, self).__init__(A=1, B=0, C=0, D=1, physicalLength=0, apertureDiameter=diameter, label=label,
+                                     frontIndex=1.0, backIndex=1.0)
 
-		Returns
-		-------
-		outputRay : object of ray class
-			the height and angle of the output ray.
+    def deviationAngle(self):
+        """ This function provides deviation angle delta assuming that
+        the axicon is in air and that the incidence is near normal,
+        which is the usual way of using an axicon.
 
-		See Also
-		--------
-		raytracing.Matrix.mul_ray
+        Returns
+        -------
+        delta : float
+            the deviation angle
 
+        See ALso
+        --------
+        https://ru.b-ok2.org/book/2482970/9062b7, p.48
+
+        """
+
+        return (self.n - 1.0) * self.alpha
+
+    def focalLineLength(self, yMax=None):
+        """ Provides the line length, assuming a ray at height yMax
+
+        Parameters
+        ----------
+        yMax : float
+            the height of the ray (default=None)
+            If no height is defined for the ray, then yMax would be set to the height of the axicon (apertureDiameter/2)
+
+        Returns
+        -------
+        focalLineLength : float
+            the length of the focal line
+
+        See ALso
+        --------
+        https://ru.b-ok2.org/book/2482970/9062b7, p.48
+
+        """
+
+        if yMax == None:
+            yMax = self.apertureDiameter / 2
+
+        return abs(yMax) / (self.n - 1.0) / self.alpha
+
+    def mul_ray(self, rightSideRay):
+        """ This function is used to calculate the output ray through an axicon.
+
+        Parameters
+        ----------
+        rightSideRay : object of ray class
+            A ray with a defined height and angle.
+
+        Returns
+        -------
+        outputRay : object of ray class
+            the height and angle of the output ray.
+
+        See Also
+        --------
+        raytracing.Matrix.mul_ray
+
+
+        """
+
+        outputRay = super(Axicon, self).mul_ray(rightSideRay)
+
+        if rightSideRay.y > 0:
+            outputRay.theta += -self.deviationAngle()
+        elif rightSideRay.y < 0:
+            outputRay.theta += self.deviationAngle()
+        # theta == 0 is not deviated
+
+        return outputRay
+
+    def mul_matrix(self, rightSideMatrix):
+        """ The final matrix of an optical path with an axicon can be calculated using this function.
+
+        Parameters
+        ----------
+        rightSideMatrix : object of matrix class
+            The ABCD matrix of an element or an optical path.
 
-		"""
+        Notes
+        -----
+        For now the final matrix with an axicon in the path cannot be calculated.
 
-		outputRay = super(Axicon, self).mul_ray(rightSideRay)
+        """
 
-		if rightSideRay.y > 0:
-			outputRay.theta += -self.deviationAngle()
-		elif rightSideRay.y < 0:
-			outputRay.theta +=  self.deviationAngle()
-		# theta == 0 is not deviated
-				
-		return outputRay
+        raise TypeError("Cannot calculate final matrix with axicon in path. \
+            You can only propagate rays all the way through")
 
-	def mul_mat(self, rightSideMatrix):
-		""" The final matrix of an optical path with an axicon can be calculated using this function.
+    def mul_beam(self, rightSideBeam):
+        """This function calculates the multiplication of a coherent beam with complex radius
+        of curvature q by an ABCD matrix.
 
-		Parameters
-		----------
-		rightSideMatrix : object of matrix class
-			The ABCD matrix of an element or an optical path.
+        Parameters
+        ----------
+        rightSideBeam : object from GaussianBeam class
+            including the beam properties
 
-		Notes
-		-----
-		For now the final matrix with an axicon in the path cannot be calculated.
 
-		"""
+        Returns
+        -------
+        outputBeam : object from GaussianBeam class
+            The properties of the beam at the output of the system with the defined ABCD matrix
 
-		raise TypeError("Cannot calculate final matrix with axicon in path. \
-			You can only propagate rays all rhe way through")
+        Examples
+        --------
+        >>> from raytracing import *
+        >>> # M1 is an ABCD matrix of a lens (f=10)
+        >>> M1= Matrix(A=1,B=0,C=-1/10,D=1,physicalLength=5,label='Lens')
+        >>> # B is a Gaussian Beam
+        >>> B=GaussianBeam(q=complex(1,1),w=1,R=5,n=1)
+        >>> print('The output properties of are:' , M1.mul_beam(B))
+        The output ray of Lens M1 : Complex radius: (0.976+1.22j)
+        w(z): 0.020, R(z): 2.500, z: 5.000, λ: 632.8 nm
+        zo: 1.220, wo: 0.016, wo position: -0.976
 
-	def drawAt(self, z, axes):
-		halfHeight = 4
-		if self.apertureDiameter != float('Inf'):
-			halfHeight = self.apertureDiameter/2
+        See Also
+        --------
+        raytracing.Matrix.mul_matrix
+        raytracing.Matrix.mul_ray
+        raytracing.GaussianBeam
+        """
 
-		plt.arrow(z, 0, 0, halfHeight, width=0.1, fc='k', ec='k',head_length=0.25, head_width=0.25,length_includes_head=True)
-		plt.arrow(z, 0, 0, -halfHeight, width=0.1, fc='k', ec='k',head_length=0.25, head_width=0.25, length_includes_head=True)
+        raise TypeError("Cannot use Axicon with GaussianBeam, only with Ray")
 
+    def drawAt(self, z, axes):  # pragma: no cover
+        halfHeight = 4
+        if self.apertureDiameter != float('Inf'):
+            halfHeight = self.apertureDiameter / 2
 
+        plt.arrow(z, 0, 0, halfHeight, width=0.1, fc='k', ec='k', head_length=0.25, head_width=0.25,
+                  length_includes_head=True)
+        plt.arrow(z, 0, 0, -halfHeight, width=0.1, fc='k', ec='k', head_length=0.25, head_width=0.25,
+                  length_includes_head=True)


Test patch:
diff --git a/raytracing/tests/testsAxicon.py b/raytracing/tests/testsAxicon.py
new file mode 100644
index 0000000..483cafa
--- /dev/null
+++ b/raytracing/tests/testsAxicon.py
@@ -0,0 +1,108 @@
+import envtest
+
+from raytracing import *
+from numpy import random
+from numpy import *
+
+inf = float("+inf")
+degrees = math.pi/180
+
+class TestAxicon(envtest.RaytracingTestCase):
+
+    def testAxicon(self):
+        n = 1.5
+        alpha = 2.6*degrees
+        diameter = 100
+        label = "Axicon"
+        axicon = Axicon(alpha, n, diameter, label)
+        self.assertEqual(axicon.n, n)
+        self.assertEqual(axicon.alpha, alpha)
+        self.assertEqual(axicon.apertureDiameter, diameter)
+        self.assertEqual(axicon.label, label)
+        self.assertEqual(axicon.frontIndex, 1.0)
+        self.assertEqual(axicon.backIndex, 1.0)
+
+    def testDeviationAngleIs0(self):
+        n = 1
+        alpha = random.randint(2000, 5000, 1).item() / 1000*degrees
+        axicon = Axicon(alpha, n)
+        self.assertEqual(axicon.deviationAngle(), 0)
+
+    def testDeviationAngleIs0Too(self):
+        n = random.randint(1000, 3000, 1).item() / 1000
+        alpha = 0*degrees
+        axicon = Axicon(alpha, n)
+        self.assertEqual(axicon.deviationAngle(), 0)
+
+    def testDeviationAngle(self):
+        n = 1.33
+        alpha = 4*degrees
+        axicon = Axicon(alpha, n)
+        self.assertAlmostEqual(axicon.deviationAngle(), 1.32*degrees, places=15)
+
+    def testFocalLineLengthYIsNoneAndInfiniteDiameter(self):
+        n = random.randint(1000, 3000, 1).item() / 1000
+        alpha = random.randint(2000, 5000, 1).item() / 1000*degrees
+        axicon = Axicon(alpha, n)
+        self.assertEqual(axicon.focalLineLength(), inf)
+
+    def testFocalLineLengthYIsNone(self):
+        n = 1.5
+        alpha = 2.6*degrees
+        axicon = Axicon(alpha, n, 100)
+        y = 50
+        L = y/tan(axicon.deviationAngle())
+        self.assertAlmostEqual(axicon.focalLineLength(), L,0)
+
+    def testFocalLineLengthSignOfY(self):
+        n = 1.43
+        alpha = 1.95*degrees
+        axicon = Axicon(alpha=alpha, n=n, diameter=100)
+        self.assertAlmostEqual(axicon.focalLineLength(-2), axicon.focalLineLength(2))
+
+    def testFocalLineLengthPositiveY(self):
+        n = 1.43
+        alpha = 1.95*degrees
+        axicon = Axicon(alpha, n, 100)
+        y = 2
+        L = y/tan(axicon.deviationAngle())
+        self.assertAlmostEqual(axicon.focalLineLength(y), L, 1)
+
+    def testHighRayIsDeviatedDown(self):
+        ray = Ray(10, 0)
+        n = 1.1
+        alpha = 2.56*degrees
+        axicon = Axicon(alpha, n, 50)
+        outputRay = axicon*ray
+        self.assertEqual(outputRay.theta, -axicon.deviationAngle())
+        self.assertTrue(outputRay.theta < 0)
+
+    def testLowRayIsDeviatedUp(self):
+        ray = Ray(-10, 0)
+        n = 1.1
+        alpha = 2.56*degrees
+        axicon = Axicon(alpha, n, 50)
+        outputRay = axicon*ray
+        self.assertEqual(outputRay.theta, axicon.deviationAngle())
+        self.assertTrue(outputRay.theta > 0)
+
+    def testMulMatrix(self):
+        matrix = Matrix()
+        axicon = Axicon(2.6543, 1.2*degrees)
+        with self.assertRaises(TypeError):
+            axicon.mul_matrix(matrix)
+
+    def testDifferentMultiplications(self):
+        ray = Ray()
+        beam = GaussianBeam(w=1, R=10, n=1.67)
+        matrix = Matrix()
+        axicon = Axicon(4.3, 1.67*degrees)
+        self.assertIsNotNone(axicon * ray)
+        with self.assertRaises(TypeError):
+            axicon * beam
+
+        with self.assertRaises(TypeError):
+            axicon * matrix
+
+if __name__ == '__main__':
+    envtest.main()
diff --git a/raytracing/tests/testsMatrixGroup.py b/raytracing/tests/testsMatrixGroup.py
index 95fa3ac..f0e835f 100644
--- a/raytracing/tests/testsMatrixGroup.py
+++ b/raytracing/tests/testsMatrixGroup.py
@@ -4,7 +4,7 @@ from raytracing import *
 
 inf = float("+inf")
 
-testSaveHugeFile = True
+testSaveHugeFiles = True
 
 
 class TestMatrixGroup(envtest.RaytracingTestCase):
@@ -591,11 +591,11 @@ class TestSaveAndLoadMatrixGroup(envtest.RaytracingTestCase):
         mg = MatrixGroup([Space(20), ThickLens(1.22, 10, 10, 10)])
         self.assertSaveNotFailed(mg, self.fileName)
 
-    @envtest.skipIf(not testSaveHugeFile, "Don't test saving a lot of matrices")
+    @envtest.skipIf(not testSaveHugeFiles, "Don't test saving a lot of matrices")
     def testSaveHugeFile(self):
         fname = self.tempFilePath("hugeFile.pkl")
-        spaces = [Space(10) for _ in range(500)]
-        lenses = [Lens(10) for _ in range(500)]
+        spaces = [Space(10) for _ in range(200)]
+        lenses = [Lens(10) for _ in range(200)]
         elements = spaces + lenses
         mg = MatrixGroup(elements)
         self.assertSaveNotFailed(mg, fname)
@@ -655,11 +655,11 @@ class TestSaveAndLoadMatrixGroup(envtest.RaytracingTestCase):
         self.assertLoadNotFailed(mg2, fname)
         self.assertLoadEqualsMatrixGroup(mg2, mg1)
 
-    @envtest.skipIf(not testSaveHugeFile, "Don't test saving a lot of matrices")
+    @envtest.skipIf(not testSaveHugeFiles, "Don't test saving a lot of matrices")
     def testSaveThenLoadHugeFile(self):
         fname = self.tempFilePath("hugeFile.pkl")
-        spaces = [Space(10) for _ in range(500)]
-        lenses = [Lens(10) for _ in range(500)]
+        spaces = [Space(10) for _ in range(125)]
+        lenses = [Lens(10) for _ in range(125)]
         elements = spaces + lenses
         mg1 = MatrixGroup(elements)
         mg2 = MatrixGroup()