Repository: DCC-Lab/RayTracing

Issue:
Matrix with incorrect determinants should not be allowed
I am working on getting everything to work properly for focal lengths when the index of refraction is not the same on either side (n1, n2) instead of 1 and 1.  Because when we do that we must make the assumption that the determinant is n1/n2, I added in the branch `differentIndicesOnEitherSides` an exception if the determinant is incorrect (which means the matrix is not physical).  This unfortunately caught many tests that were not really correct, and now (in that branch) we get 42 errors and 6 failures when running all tests. I need help fixing those tests. @GabGabG can you help?

Discussion:

Patch:
diff --git a/.idea/other.xml b/.idea/other.xml
new file mode 100644
index 0000000..640fd80
--- /dev/null
+++ b/.idea/other.xml
@@ -0,0 +1,7 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="PySciProjectComponent">
+    <option name="PY_SCI_VIEW" value="true" />
+    <option name="PY_SCI_VIEW_SUGGESTED" value="true" />
+  </component>
+</project>
\ No newline at end of file
diff --git a/raytracing/matrix.py b/raytracing/matrix.py
index 10b2044..2f8ff11 100644
--- a/raytracing/matrix.py
+++ b/raytracing/matrix.py
@@ -135,6 +135,8 @@ class Matrix(object):
         self.label = label
         self.isFlipped = False
         super(Matrix, self).__init__()
+        if not isclose(self.determinant, self.frontIndex / self.backIndex, atol=self.__epsilon__):
+            raise ValueError("The matrix has inconsistent values")
 
     @property
     def isIdentity(self):
@@ -143,7 +145,10 @@ class Matrix(object):
     @property
     def determinant(self):
         """The determinant of the ABCD matrix is always frontIndex/backIndex,
-        which is often 1.0
+        which is often 1.0. 
+        We make a calculation exception when C == 0 and B is infinity: since
+        B is never really infinity, but C can be precisely zero (especially
+        in free space), then B*C is zero in that particular case.
 
         Examples
         --------
@@ -155,6 +160,10 @@ class Matrix(object):
         the determinant of matrix is equal to : 1.0
 
         """
+
+        if self.C == 0:
+            return self.A * self.D
+        
         return self.A * self.D - self.B * self.C
 
     def __mul__(self, rightSide):
@@ -921,10 +930,18 @@ class Matrix(object):
         """ The effective focal lengths calculated from the power (C)
         of the matrix.
 
+        There are in general two effective focal lengths: front effective
+        and back effective focal lengths (not to be confused with back focal
+        and front focal lengths which are measured from the physical interface).
+        The easiest way to calculate this is to use
+        f = -1/C for current matrix, then flipOrientation and f = -1/C
+
+
         Returns
         -------
         effectiveFocalLengths : array
-            Returns the FFL and BFL
+            Returns the effective focal lengths in the forward and backward
+            directions. When in air, both are equal.
 
         See Also
         --------
@@ -944,17 +961,15 @@ class Matrix(object):
         >>> print('focal distances:' , f1)
         focal distances: (5.0, 5.0)
 
-        Notes
-        -----
-        Currently, it is assumed the index is n=1 on either side and
-        both focal lengths are the same.
         """
         if self.hasPower:
-            focalLength = -1.0 / self.C  # FIXME: Assumes n=1 on either side
+            focalLength2 = -1.0 / self.C  # left (n1) to right (n2)
+            focalLength1 = -(self.frontIndex / self.backIndex) / self.C  # right (n2) to left (n2)
         else:
-            focalLength = float("+Inf")
+            focalLength1 = float("+Inf")
+            focalLength2 = float("+Inf")
 
-        return (focalLength, focalLength)
+        return (focalLength1, focalLength2)
 
     def backFocalLength(self):
         """ The focal lengths measured from the back vertex.
@@ -995,8 +1010,8 @@ class Matrix(object):
         we may not know where the front and back vertices are. In that case,
         we return None (or undefined).
 
-        Currently, it is assumed the index is n=1 on either side and
-        both focal distances are the same.
+        The front and back focal lengths will be different if the index
+        of refraction is different on both sides.
         """
 
         if self.backVertex is not None and self.hasPower:
@@ -1047,8 +1062,8 @@ class Matrix(object):
         we may not know where the front and back vertices are. In that case,
         we return None (or undefined).
 
-        Currently, it is assumed the index is n=1 on either side and
-        both focal distances are the same.
+        The front and back focal lengths will be different if the index
+        of refraction is different on both sides.
         """
 
         if self.frontVertex is not None and self.hasPower:
@@ -1062,10 +1077,13 @@ class Matrix(object):
     def focusPositions(self, z):
         """ Positions of both focal points on either side of the element.
 
+        The front and back focal spots will be different if the index
+        of refraction is different on both sides.
+
         Parameters
         ----------
         z : float
-            The position in which the object is placed
+            Position from where the positions are calculated
 
         Returns
         -------
@@ -1106,7 +1124,7 @@ class Matrix(object):
         Parameters
         ----------
         z : float
-            The position
+            Position from where the positions are calculated
 
         Returns
         -------
@@ -1128,8 +1146,7 @@ class Matrix(object):
         raytracing.Matrix.focusPositions
         """
         if self.hasPower:
-            p1 = z - (1 - self.D) / self.C  # FIXME: Assumes n=1 on either side
-            # FIXME: Assumes n=1 on either side
+            p1 = z - (self.frontIndex / self.backIndex - self.D) / self.C
             p2 = z + self.L + (1 - self.A) / self.C
         else:
             p1 = None
@@ -1186,7 +1203,7 @@ class Matrix(object):
             conjugateMatrix = None  # Unable to compute with inf
         else:
             distance = -self.B / self.D
-            conjugateMatrix = Space(d=distance) * self
+            conjugateMatrix = Space(d=distance, n=self.backIndex) * self
 
         return (distance, conjugateMatrix)
 


Test patch:
diff --git a/raytracing/tests/testsMatrix.py b/raytracing/tests/testsMatrix.py
index be36c0a..f36fe01 100644
--- a/raytracing/tests/testsMatrix.py
+++ b/raytracing/tests/testsMatrix.py
@@ -1,5 +1,7 @@
 import envtest  # modifies path
+import platform
 import sys
+
 from raytracing import *
 
 inf = float("+inf")
@@ -19,22 +21,26 @@ class TestMatrix(envtest.RaytracingTestCase):
         self.assertIsNotNone(m)
 
     def testMatrixExplicit(self):
-        m = Matrix(A=1, B=2, C=3, D=4, physicalLength=1,
+        m = Matrix(A=1, B=0, C=0, D=1, physicalLength=1,
                    frontVertex=0, backVertex=0, apertureDiameter=1.0)
         self.assertIsNotNone(m)
         self.assertEqual(m.A, 1)
-        self.assertEqual(m.B, 2)
-        self.assertEqual(m.C, 3)
-        self.assertEqual(m.D, 4)
+        self.assertEqual(m.B, 0)
+        self.assertEqual(m.C, 0)
+        self.assertEqual(m.D, 1)
+        self.assertEqual(m.L, 1)
+        self.assertEqual(m.backVertex, 0)
+        self.assertEqual(m.frontVertex, 0)
+        self.assertEqual(m.apertureDiameter, 1)
 
     def testMatrixProductMath(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4)
-        m2 = Matrix(A=5, B=6, C=7, D=8)
+        m1 = Matrix(A=4, B=3, C=1, D=1)
+        m2 = Matrix(A=1, B=1, C=3, D=4)
         m3 = m2 * m1
-        self.assertEqual(m3.A, 1 * 5 + 3 * 6)
-        self.assertEqual(m3.B, 2 * 5 + 4 * 6)
-        self.assertEqual(m3.C, 1 * 7 + 3 * 8)
-        self.assertEqual(m3.D, 2 * 7 + 4 * 8)
+        self.assertEqual(m3.A, 5)
+        self.assertEqual(m3.B, 4)
+        self.assertEqual(m3.C, 16)
+        self.assertEqual(m3.D, 13)
 
     def testIsIdentity(self):
         m = Matrix()
@@ -52,117 +58,117 @@ class TestMatrix(envtest.RaytracingTestCase):
         self.assertEqual(m3.backIndex, 1)
 
     def testMatrixProductIndicesLHSIsIdentity(self):
-        m1 = Matrix(backIndex=1.33)
+        m1 = Matrix()
         m2 = Matrix(1, 10, 0, 1, frontIndex=1.5, backIndex=1.5)
         m3 = m1 * m2
         self.assertEqual(m3.frontIndex, 1.5)
         self.assertEqual(m3.backIndex, 1.5)
 
     def testMatrixProductIndicesRHSIsIdentity(self):
-        m1 = Matrix(backIndex=1.33)
+        m1 = Matrix()
         m2 = Matrix(1, 10, 0, 1, frontIndex=1.5, backIndex=1.5)
         m3 = m2 * m1
         self.assertEqual(m3.frontIndex, 1.5)
         self.assertEqual(m3.backIndex, 1.5)
 
     def testMatrixProductIndicesNoIdentity(self):
-        m1 = Matrix(1, 10, 0, 1, backIndex=1.33, frontIndex=1)
-        m2 = Matrix(1, 10, 0, 1, backIndex=1, frontIndex=1.33)
+        m1 = Matrix(1, 10, 0, 0.7518796992, backIndex=1.33, frontIndex=1)
+        m2 = Matrix(1.33, 10, 0, 1, backIndex=1, frontIndex=1.33)
         m3 = m2 * m1
         self.assertEqual(m3.frontIndex, 1)
         self.assertEqual(m3.backIndex, 1)
 
     def testMatrixProductWithRayMath(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4)
+        m1 = Matrix(A=1, B=1, C=3, D=4)
         rayIn = Ray(y=1, theta=0.1)
         rayOut = m1 * rayIn
-        self.assertEqual(rayOut.y, 1 * 1 + 2 * 0.1)
-        self.assertEqual(rayOut.theta, 3 * 1 + 4 * 0.1)
+        self.assertEqual(rayOut.y, 1.1)
+        self.assertEqual(rayOut.theta, 3.4)
 
-    def testMatrixProductOutpuRayLength(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, physicalLength=2)
+    def testMatrixProductOutputRayLength(self):
+        m1 = Matrix(A=1, B=0, C=0, D=1, physicalLength=2)
         rayIn = Ray(y=1, theta=0.1, z=1)
         rayOut = m1 * rayIn
         self.assertEqual(rayOut.z, 2 + 1)
 
     def testMatrixProductOutputRayAperture(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, physicalLength=2)
+        m1 = Matrix(A=1, B=0, C=0, D=1, physicalLength=2)
         rayIn = Ray(y=1, theta=0.1, z=1)
         rayOut = m1 * rayIn
         self.assertEqual(rayOut.apertureDiameter, inf)
 
     def testMatrixProductWithRayGoesOverAperture(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, apertureDiameter=10)
+        m1 = Matrix(A=1, B=0, C=0, D=1, apertureDiameter=10)
         rayIn = Ray(y=6, theta=0.1, z=1)
         rayOut = m1 * rayIn
         self.assertTrue(rayOut.isBlocked)
 
     def testMatrixProductWithRayGoesUnderAperture(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, apertureDiameter=10)
+        m1 = Matrix(A=1, B=0, C=0, D=1, apertureDiameter=10)
         rayIn = Ray(y=-6, theta=0.1, z=1)
         rayOut = m1 * rayIn
         self.assertTrue(rayOut.isBlocked)
 
     def testMatrixProductRayGoesInAperture(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, apertureDiameter=10)
+        m1 = Matrix(A=1, B=0, C=0, D=1, apertureDiameter=10)
         rayIn = Ray(y=-1, theta=0.1, z=1)
         rayOut = m1 * rayIn
         self.assertFalse(rayOut.isBlocked)
 
     def testMatrixProductRayAlreadyBlocked(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, apertureDiameter=10)
+        m1 = Matrix(A=1, B=0, C=0, D=1, apertureDiameter=10)
         rayIn = Ray(y=-1, theta=0.1, z=1, isBlocked=True)
         rayOut = m1 * rayIn
         self.assertTrue(rayOut.isBlocked)
 
     def testMatrixProductLength(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4)
-        m2 = Matrix(A=5, B=6, C=7, D=8)
+        m1 = Matrix(A=1, B=0, C=0, D=1)
+        m2 = Matrix(A=1, B=0, C=0, D=1)
         m3 = m2 * m1
         self.assertEqual(m3.L, m1.L + m2.L)
         self.assertIsNone(m3.frontVertex)
         self.assertIsNone(m3.backVertex)
 
     def testMatrixProductVertices(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, physicalLength=10, frontVertex=0, backVertex=10)
+        m1 = Matrix(A=1, B=0, C=0, D=1, physicalLength=10, frontVertex=0, backVertex=10)
         self.assertEqual(m1.frontVertex, 0)
         self.assertEqual(m1.backVertex, 10)
 
     def testMatrixProductVerticesAllNone(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4)
-        m2 = Matrix(A=5, B=6, C=7, D=8)
+        m1 = Matrix(A=1, B=0, C=0, D=1)
+        m2 = Matrix(A=1, B=0, C=0, D=1)
         m3 = m2 * m1
         self.assertEqual(m3.L, m1.L + m2.L)
         self.assertIsNone(m3.frontVertex)
         self.assertIsNone(m3.backVertex)
 
     def testMatrixProductVerticesSecondNone(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, physicalLength=10, frontVertex=0, backVertex=10)
-        m2 = Matrix(A=5, B=6, C=7, D=8)
+        m1 = Matrix(A=1, B=0, C=0, D=1, physicalLength=10, frontVertex=0, backVertex=10)
+        m2 = Matrix(A=1, B=0, C=0, D=1)
         m3 = m2 * m1
         self.assertEqual(m3.L, m1.L + m2.L)
         self.assertEqual(m3.frontVertex, 0)
         self.assertEqual(m3.backVertex, 10)
 
     def testMatrixProductVerticesFirstNone(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4)
-        m2 = Matrix(A=5, B=6, C=7, D=8, physicalLength=10, frontVertex=0, backVertex=10)
+        m1 = Matrix(A=1, B=0, C=0, D=1)
+        m2 = Matrix(A=1, B=0, C=0, D=1, physicalLength=10, frontVertex=0, backVertex=10)
         m3 = m2 * m1
         self.assertEqual(m3.L, m1.L + m2.L)
         self.assertEqual(m3.frontVertex, 0)
         self.assertEqual(m3.backVertex, 10)
 
     def testMatrixProductVerticesTwoElements(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, physicalLength=5, frontVertex=0, backVertex=5)
-        m2 = Matrix(A=5, B=6, C=7, D=8, physicalLength=10, frontVertex=0, backVertex=10)
+        m1 = Matrix(A=1, B=0, C=0, D=1, physicalLength=5, frontVertex=0, backVertex=5)
+        m2 = Matrix(A=1, B=0, C=0, D=1, physicalLength=10, frontVertex=0, backVertex=10)
         m3 = m2 * m1
         self.assertEqual(m3.L, m1.L + m2.L)
         self.assertEqual(m3.frontVertex, 0)
         self.assertEqual(m3.backVertex, 15)
 
     def testMatrixProductVerticesTwoElementsRepresentingGroups(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, physicalLength=5, frontVertex=1, backVertex=4)
-        m2 = Matrix(A=5, B=6, C=7, D=8, physicalLength=10, frontVertex=2, backVertex=9)
+        m1 = Matrix(A=1, B=0, C=0, D=1, physicalLength=5, frontVertex=1, backVertex=4)
+        m2 = Matrix(A=1, B=0, C=0, D=1, physicalLength=10, frontVertex=2, backVertex=9)
         m3 = m2 * m1
         self.assertEqual(m3.L, m1.L + m2.L)
         self.assertEqual(m3.frontVertex, 1)
@@ -174,14 +180,13 @@ class TestMatrix(envtest.RaytracingTestCase):
         self.assertEqual(m3.backVertex, 14)
 
     def testMatrixProductGaussianBeamMath(self):
-        m = Matrix(A=1, B=2, C=3, D=4)
-        beamIn = GaussianBeam(w=1, wavelength=1)  # q = j\pi
+        m = Matrix(A=2, B=1, C=3, D=2)
+        beamIn = GaussianBeam(w=1, wavelength=1)  # q = j*pi
         beamOut = m * beamIn
-        q = complex(0, math.pi)
-        self.assertEqual(beamOut.q, (1 * q + 2) / (3 * q + 4))
+        self.assertEqual(beamOut.q, (2j * pi + 1) / (3j * pi + 2))
 
     def testMatrixProductGaussianNotSameRefractionIndex(self):
-        m = Matrix(A=1, B=2, C=3, D=4)
+        m = Matrix(A=1, B=0, C=0, D=1)
         beam = GaussianBeam(w=1, n=1.2)
 
         with self.assertRaises(UserWarning):
@@ -190,38 +195,38 @@ class TestMatrix(envtest.RaytracingTestCase):
                 m * beam
 
     def testMatrixProductGaussianBeamWavelengthOut(self):
-        m = Matrix(A=1, B=2, C=3, D=4, )
+        m = Matrix(A=1, B=0, C=0, D=1, )
         beamIn = GaussianBeam(w=1, wavelength=1)
         beamOut = m * beamIn
         self.assertEqual(beamOut.wavelength, 1)
 
     def testMatrixProductGaussianRefractIndexOut(self):
-        m = Matrix(A=1, B=2, C=3, D=4, frontIndex=1.33, backIndex=1.33)
+        m = Matrix(A=1, B=0, C=0, D=1, frontIndex=1.33, backIndex=1.33)
         beamIn = GaussianBeam(w=1, wavelength=1, n=1.33)
         beamOut = m * beamIn
         self.assertEqual(beamOut.n, 1.33)
 
     def testMatrixProductGaussianLength(self):
-        m = Matrix(A=1, B=2, C=3, D=4, frontIndex=1.33, physicalLength=1.2)
+        m = Matrix(A=1, B=0, C=0, D=1.33, frontIndex=1.33, physicalLength=1.2)
         beamIn = GaussianBeam(w=1, wavelength=1, z=1, n=1.33)
         beamOut = m * beamIn
         self.assertEqual(beamOut.z, 2.2)
 
     def testMatrixProductGaussianClippedOverAperture(self):
-        m = Matrix(A=1, B=2, C=3, D=4, physicalLength=1.2, apertureDiameter=2)
+        m = Matrix(A=1, B=0, C=0, D=1, physicalLength=1.2, apertureDiameter=2)
         beamIn = GaussianBeam(w=1.1, wavelength=1, z=1)
         beamOut = m * beamIn
         self.assertTrue(beamOut.isClipped)
 
     def testMatrixProductGaussianInitiallyClipped(self):
-        m = Matrix(A=1, B=2, C=3, D=4, physicalLength=1.2, apertureDiameter=2)
+        m = Matrix(A=1, B=0, C=0, D=1, physicalLength=1.2, apertureDiameter=2)
         beamIn = GaussianBeam(w=0.5, wavelength=1, z=1)
         beamIn.isClipped = True
         beamOut = m * beamIn
         self.assertTrue(beamOut.isClipped)
 
     def testMatrixProductGaussianNotClipped(self):
-        m = Matrix(A=1, B=2, C=3, D=4, physicalLength=1.2)
+        m = Matrix(A=1, B=0, C=0, D=1, physicalLength=1.2)
         beamIn = GaussianBeam(w=1.1, wavelength=1, z=1)
         beamOut = m * beamIn
         self.assertFalse(beamOut.isClipped)
@@ -233,20 +238,20 @@ class TestMatrix(envtest.RaytracingTestCase):
             m * other
 
     def testApertureDiameter(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, apertureDiameter=2)
+        m1 = Matrix(A=1, B=0, C=0, D=1, apertureDiameter=2)
         self.assertTrue(m1.hasFiniteApertureDiameter())
         self.assertEqual(m1.largestDiameter, 2.0)
-        m2 = Matrix(A=1, B=2, C=3, D=4)
+        m2 = Matrix(A=1, B=0, C=0, D=1)
         self.assertFalse(m2.hasFiniteApertureDiameter())
         self.assertEqual(m2.largestDiameter, float("+inf"))
 
     def testTransferMatrix(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4)
+        m1 = Matrix(A=1, B=0, C=0, D=1)
         # Null length returns self
         self.assertEqual(m1.transferMatrix(), m1)
 
         # Length == 1 returns self if upTo >= 1
-        m2 = Matrix(A=1, B=2, C=3, D=4, physicalLength=1)
+        m2 = Matrix(A=1, B=0, C=0, D=1, physicalLength=1)
         self.assertEqual(m2.transferMatrix(upTo=1), m2)
         self.assertEqual(m2.transferMatrix(upTo=2), m2)
 
@@ -257,31 +262,30 @@ class TestMatrix(envtest.RaytracingTestCase):
             m2.transferMatrix(upTo=0.5)
 
     def testTransferMatrices(self):
-        m1 = Matrix(A=1, B=2, C=3, D=4, frontIndex=2)
+        m1 = Matrix(A=1, B=0, C=0, D=2, frontIndex=2)
         self.assertEqual(m1.transferMatrices(), [m1])
-        m1 * GaussianBeam(w=1, n=2)
 
     def testTrace(self):
         ray = Ray(y=1, theta=1)
-        m = Matrix(A=1, B=2, C=3, D=4, physicalLength=1)
+        m = Matrix(A=1, B=0, C=0, D=1, physicalLength=1)
         trace = [ray, m * ray]
         self.assertListEqual(m.trace(ray), trace)
 
     def testTraceNullLength(self):
         ray = Ray(y=1, theta=1)
-        m = Matrix(A=1, B=2, C=3, D=4)
+        m = Matrix(A=1, B=0, C=0, D=1)
         trace = [m * ray]
         self.assertListEqual(m.trace(ray), trace)
 
     def testTraceBlocked(self):
         ray = Ray(y=10, theta=1)
-        m = Matrix(A=1, B=2, C=3, D=4, apertureDiameter=10, physicalLength=1)
+        m = Matrix(A=1, B=0, C=0, D=1, apertureDiameter=10, physicalLength=1)
         trace = m.trace(ray)
         self.assertTrue(all(x.isBlocked for x in trace))
 
     def testTraceGaussianBeam(self):
         beam = GaussianBeam(w=1)
-        m = Matrix(A=1, B=2, C=3, D=4, apertureDiameter=10)
+        m = Matrix(A=1, B=0, C=0, D=1, apertureDiameter=10)
         outputBeam = m * beam
         tracedBeam = m.trace(beam)[-1]
         self.assertEqual(tracedBeam.w, outputBeam.w)
@@ -292,7 +296,7 @@ class TestMatrix(envtest.RaytracingTestCase):
 
     def testTraceThrough(self):
         ray = Ray()
-        m = Matrix(A=1, B=2, C=3, D=4, apertureDiameter=10)
+        m = Matrix(A=1, B=0, C=0, D=1, apertureDiameter=10)
         trace = m.traceThrough(ray)
         self.assertEqual(trace, m * ray)
 
@@ -383,16 +387,16 @@ class TestMatrix(envtest.RaytracingTestCase):
             traceWithNumberProcesses = m.traceManyThroughInParallel(rays, processes=2)
             for i in range(len(rays)):
                 # Order is not kept, we have to check if the ray traced is in the original list
-                self.assertTrue(traceWithNumberProcesses[i] in rays)
-        except:
-            pass
+                self.assertIn(traceWithNumberProcesses[i], rays)
+        except Exception as exception:
+            self.fail(f"Exception raised:\n{exception}")
 
     def testPointsOfInterest(self):
         m = Matrix()
         self.assertListEqual(m.pointsOfInterest(1), [])
 
     def testIsImaging(self):
-        m = Matrix(A=1, B=0, C=3, D=4)
+        m = Matrix(A=1, B=0, C=3, D=1)
         self.assertTrue(m.isImaging)
 
     def testIsNotImaging(self):
@@ -416,39 +420,41 @@ class TestMatrix(envtest.RaytracingTestCase):
         self.assertFalse(m.hasPower)
 
     def testEffectiveFocalLengthsHasPower(self):
-        m = Matrix(1, 2, 3, 4)
+        m = Matrix(A=1, B=0, C=3, D=1)
         focalLengths = (-1 / 3, -1 / 3)
         self.assertTupleEqual(m.effectiveFocalLengths(), focalLengths)
 
     def testEffectiveFocalLengthsNoPower(self):
-        m = Matrix()
+        m = Matrix(1, 0, 0, 1)
         focalLengths = (inf, inf)
         self.assertTupleEqual(m.effectiveFocalLengths(), focalLengths)
 
     def testMatrixBackFocalLength(self):
-        m = Matrix(1, 2, 3, 4, backVertex=1, physicalLength=1)
-        f2 = -1 / 3
-        p2 = 0 + 1 + (1 - 1) / 3
-        self.assertEqual(m.backFocalLength(), p2 + f2 - 1)
+        R = 10
+        n1 = 1.2
+        n2 = 1.5
+        m = DielectricInterface(n1=n1, n2=n2, R=R)
+        self.assertAlmostEqual(m.backFocalLength(), -m.A / m.C)
 
     def testBackFocalLengthSupposedNone(self):
         m = Matrix()
         self.assertIsNone(m.backFocalLength())
 
     def testMatrixFrontFocalLength(self):
-        m = Matrix(1, 2, 3, 4, frontVertex=1, physicalLength=1)
-        f1 = -1 / 3
-        p1 = 0 - (1 - 4) / 3
-        self.assertEqual(m.frontFocalLength(), -(p1 - f1 - 1))
+        R = 10
+        n1 = 1.2
+        n2 = 1.5
+        m = DielectricInterface(n1=n1, n2=n2, R=R)
+        self.assertAlmostEqual(m.frontFocalLength(), -m.D / m.C)
 
     def testFrontFocalLengthSupposedNone(self):
         m = Matrix()
         self.assertIsNone(m.frontFocalLength())
 
     def testPrincipalPlanePositions(self):
-        m = Matrix(1, 2, 3, 4, physicalLength=1)
-        p1 = 0 - (1 - 4) / 3
-        p2 = 0 + 1 + (1 - 1) / 3
+        m = Matrix(A=1, B=0, C=1, D=1, physicalLength=1)
+        p1 = 0
+        p2 = 1
         self.assertTupleEqual(m.principalPlanePositions(0), (p1, p2))
 
     def testPrincipalPlanePositionsNoPower(self):
@@ -456,19 +462,29 @@ class TestMatrix(envtest.RaytracingTestCase):
         self.assertTupleEqual(m.principalPlanePositions(0), (None, None))
 
     def testFocusPositions(self):
-        m = Matrix(1, 2, 3, 4, physicalLength=1)
-        f1 = -1 / 3
-        p1 = 1
-        f2 = -1 / 3
-        p2 = 1
+        m = Matrix(A=1 / 3, B=0, C=10, D=3, physicalLength=1)
+        f1 = -0.1
+        p1 = 0.2
+        f2 = -0.1
+        p2 = 16 / 15
         self.assertTupleEqual(m.focusPositions(0), (p1 - f1, p2 + f2))
 
     def testFocusPositionsNoPower(self):
         m = Matrix()
         self.assertTupleEqual(m.focusPositions(0), (None, None))
 
-    def testFiniteForwardConjugate_1(self):
-        m1 = Matrix(1, 0, -1 / 5, 1) * Matrix(1, 10, 0, 1)
+    def testDielectricInterfaceEffectiveFocalLengths(self):
+        # Positive R is convex for ray
+        n1 = 1
+        n2 = 1.5
+        R = 10
+        m = DielectricInterface(n1=n1, n2=n2, R=R)
+        (f1, f2) = m.effectiveFocalLengths()
+        self.assertTrue(f2 == n2 * R / (n2 - n1))
+        self.assertTrue(f1 == n1 * R / (n2 - n1))  # flip R and n1,n2
+
+    def testFiniteForwardConjugate(self):
+        m1 = Lens(f=5) * Space(d=10)
         (d, m2) = m1.forwardConjugate()
         self.assertTrue(m2.isImaging)
         self.assertEqual(d, 10)
@@ -486,12 +502,12 @@ class TestMatrix(envtest.RaytracingTestCase):
         m1 = Matrix(1, 0, -1 / 5, 1) * Matrix(1, 5, 0, 1)
         (d, m2) = m1.forwardConjugate()
         self.assertIsNone(m2)
-        self.assertEqual(d, float("+inf"))
+        self.assertEqual(d, inf)
         self.assertEqual(m1.determinant, 1)
 
     def testInfiniteBackConjugate(self):
-        m = Matrix(A=0)
-        self.assertTupleEqual(m.backwardConjugate(), (float("+inf"), None))
+        m = Matrix(A=0, B=1, C=-1)
+        self.assertTupleEqual(m.backwardConjugate(), (inf, None))
 
     def testFiniteBackConjugate_1(self):
         m1 = Matrix(1, 10, 0, 1) * Matrix(1, 0, -1 / 5, 1)
@@ -522,7 +538,7 @@ class TestMatrix(envtest.RaytracingTestCase):
         backVertexInit = 20
         frontIndexInit = 1
         backIndexInit = 2
-        m = Matrix(frontVertex=frontVertexInit, backVertex=backVertexInit, frontIndex=frontIndexInit,
+        m = Matrix(A=0.5, frontVertex=frontVertexInit, backVertex=backVertexInit, frontIndex=frontIndexInit,
                    backIndex=backIndexInit)
         m.flipOrientation()
         self.assertTrue(m.isFlipped)
@@ -570,12 +586,16 @@ class TestMatrix(envtest.RaytracingTestCase):
         self.assertNotEqual(m, "Trust me, this is a Matrix. This is equal to Matrix()")
 
     def testEqualityMatricesNotEqualSameABCD(self):
-        m = Matrix()
-        m2 = Matrix(frontIndex=10)
+        m = Matrix(1,0,0,1)
+        m2 = Matrix(1,0,0,1, frontVertex=1)
+        self.assertNotEqual(m, m2)
+        m2 = Matrix(1,0,0,1, backVertex=1)
+        self.assertNotEqual(m, m2)
+        m2 = Matrix(1,0,0,1, frontIndex=10, backIndex=10)
         self.assertNotEqual(m, m2)
 
     def testEqualityMatricesNotEqualDifferentABCD(self):
-        m = Matrix()
+        m = Matrix(1,0,0,1)
         m2 = Matrix(A=1 / 2, D=2)
         self.assertNotEqual(m, m2)
 
diff --git a/raytracing/tests/testsMatrixGroup.py b/raytracing/tests/testsMatrixGroup.py
index 4d843ce..95fa3ac 100644
--- a/raytracing/tests/testsMatrixGroup.py
+++ b/raytracing/tests/testsMatrixGroup.py
@@ -252,7 +252,7 @@ class TestMatrixGroup(envtest.RaytracingTestCase):
         self.assertListEqual(mg.intermediateConjugates(), [])
 
     def testIntermediateConjugatesNoConjugate(self):
-        mg = MatrixGroup([Matrix(1, 1, 1, 0, 1)])
+        mg = MatrixGroup([Matrix(1, 1, -1, 0, 1)])
         self.assertListEqual(mg.intermediateConjugates(), [])
 
     def testIntermediateConjugates(self):