Repository: DCC-Lab/RayTracing

Issue:
demo #3-13-16-18-19 not working with 1.2.10
I installed raytracing with pip.
i called "python -m raytracing"
demo 1 and 2 worked, but not 3, so others didnt display.
i called every demo 1 by 1 to see which one is working.
It looks like some demos are not displaying correclty (i think).



Discussion:

Patch:
diff --git a/raytracing/__main__.py b/raytracing/__main__.py
index 1db0a3d..dd34f4f 100644
--- a/raytracing/__main__.py
+++ b/raytracing/__main__.py
@@ -1,5 +1,6 @@
 from .imagingpath import *
 from .laserpath import *
+from .lasercavity import *
 
 from .specialtylenses import *
 from .axicon import *
@@ -392,7 +393,7 @@ if 18 in examples:
     path.append(Space(d=180))
     path.append(olympus.LUMPlanFL40X())
     path.append(Space(d=10))
-    path.display(inputBeam=GaussianBeam(w=0.001), comments="""
+    path.display(beams=[GaussianBeam(w=0.001)], comments="""
     path = LaserPath()
     path.label = "Demo #18: Laser beam and vendor lenses"
     path.append(Space(d=50))
@@ -408,21 +409,20 @@ if 18 in examples:
     path.append(Space(d=10))
     path.display()""")
 if 19 in examples:
-    cavity = LaserPath(label="Laser cavity: round trip\nCalculated laser modes")
-    cavity.isResonator = True
+    cavity = LaserCavity(label="Laser cavity: round trip\nCalculated laser modes")
     cavity.append(Space(d=160))
     cavity.append(DielectricSlab(thickness=100, n=1.8))
     cavity.append(Space(d=160))
-    cavity.append(CurvedMirror(R=400))
+    cavity.append(CurvedMirror(R=-400))
     cavity.append(Space(d=160))
     cavity.append(DielectricSlab(thickness=100, n=1.8))
     cavity.append(Space(d=160))
 
     # Calculate all self-replicating modes (i.e. eigenmodes)
     (q1, q2) = cavity.eigenModes()
-    print(q1, q2)
+    print(q1,q2)
 
-    # Obtain all physical (i.e. finite) self-replicating modes
+    # Obtain all physical modes (i.e. only finite eigenmodes)
     qs = cavity.laserModes()
     for q in qs:
         print(q)
diff --git a/raytracing/figure.py b/raytracing/figure.py
index 6acd631..bea4212 100644
--- a/raytracing/figure.py
+++ b/raytracing/figure.py
@@ -14,6 +14,7 @@ class Figure:
         self.figure = None
         self.axes = None  # Where the optical system is
         self.axesComments = None  # Where the comments are (for teaching)
+        self.elementGraphics = []
 
         self.styles = dict()
         self.styles['default'] = {'rayColors': ['b', 'r', 'g'], 'onlyAxialRay': False,
@@ -64,7 +65,7 @@ class Figure:
 
         if self.designParams['onlyPrincipalAndAxialRays']:
             (stopPosition, stopDiameter) = self.path.apertureStop()
-            if stopPosition is None:
+            if stopPosition is None or self.path.principalRay() is None:
                 warnings.warn("No aperture stop in system: cannot use onlyPrincipalAndAxialRays=True since they are "
                               "not defined.")
                 self.designParams['onlyPrincipalAndAxialRays'] = False
@@ -136,14 +137,15 @@ class Figure:
         self.drawDisplayObjects()
 
         self.axes.callbacks.connect('ylim_changed', self.onZoomCallback)
-        self.axes.set_ylim([-self.displayRange() / 2 * 1.5, self.displayRange() / 2 * 1.5])
+        self.axes.set_xlim(0 - self.path.L * 0.05, self.path.L + self.path.L * 0.05)
+        self.axes.set_ylim([-self.displayRange() / 2 * 1.6, self.displayRange() / 2 * 1.6])
 
         if filepath is not None:
             self.figure.savefig(filepath, dpi=600)
         else:
             self._showPlot()
 
-    def displayGaussianBeam(self, inputBeams=None, filepath=None):
+    def displayGaussianBeam(self, beams=None, filepath=None):
         """ Display the optical system and trace the laser beam.
         If comments are included they will be displayed on a
         graph in the bottom half of the plot.
@@ -154,11 +156,14 @@ class Figure:
             A list of Gaussian beams
         """
 
-        self.drawBeamTraces(beams=inputBeams)
+        if len(beams) != 0:
+            self.drawBeamTraces(beams=beams)
+
         self.drawDisplayObjects()
 
         self.axes.callbacks.connect('ylim_changed', self.onZoomCallback)
-        self.axes.set_ylim([-self.displayRange() / 2 * 1.5, self.displayRange() / 2 * 1.5])
+        self.axes.set_xlim(0 - self.path.L * 0.05, self.path.L + self.path.L * 0.05)
+        self.axes.set_ylim([-self.displayRange() / 2 * 1.6, self.displayRange() / 2 * 1.6])
 
         if filepath is not None:
             self.figure.savefig(filepath, dpi=600)
@@ -272,10 +277,13 @@ class Figure:
             return self.imagingDisplayRange()
 
     def imagingDisplayRange(self):
-        displayRange = self.path.largestDiameter
+        displayRange = 0
+        for graphic in self.elementGraphics:
+            if graphic.halfHeight * 2 > displayRange:
+                displayRange = graphic.halfHeight * 2
 
-        if displayRange == float('+Inf') or displayRange <= 2 * self.path._objectHeight:
-            displayRange = 2 * self.path._objectHeight
+        if displayRange == float('+Inf') or displayRange <= self.path._objectHeight:
+            displayRange = self.path._objectHeight
 
         conjugates = self.path.intermediateConjugates()
         if len(conjugates) != 0:
@@ -289,9 +297,16 @@ class Figure:
         return displayRange
 
     def laserDisplayRange(self):
-        displayRange = self.path.largestDiameter
-        if displayRange == float('+Inf'):
-            displayRange = self.path.inputBeam.w * 3
+        displayRange = 0
+        for graphic in self.elementGraphics:
+            if graphic.halfHeight * 2 > displayRange:
+                displayRange = graphic.halfHeight * 2
+
+        if displayRange == float('+Inf') or displayRange == 0:
+            if self.path.inputBeam is not None:
+                displayRange = self.path.inputBeam.w * 3
+            else:
+                displayRange = 100
 
         return displayRange
 
@@ -515,6 +530,7 @@ class Figure:
                 color='r'))
 
     def drawElements(self, elements):
+        self.elementGraphics = []
         z = 0
         for element in elements:
             graphic = Graphic(element)
@@ -524,6 +540,7 @@ class Figure:
             if self.path.showElementLabels:
                 graphic.drawLabels(z, self.axes)
             z += graphic.L
+            self.elementGraphics.append(graphic)
 
     def rayTraceLines(self, removeBlockedRaysCompletely=True):
         """ A list of all ray trace line objects corresponding to either
@@ -613,8 +630,8 @@ class Figure:
         yScale : float
             The scale of y axes
         """
-
-        xScale, yScale = self.axes.viewLim.bounds[2:]
+        xScale = self.path.L * 1.1
+        yScale = self.displayRange() * 1.6
 
         return xScale, yScale
 
@@ -638,11 +655,18 @@ class Figure:
 class MatrixGraphic:
     def __init__(self, matrix: Matrix):
         self.matrix = matrix
+        self._halfHeight = None
 
     @property
     def L(self):
         return self.matrix.L
 
+    @property
+    def halfHeight(self):
+        if self._halfHeight is None:
+            self._halfHeight = self.displayHalfHeight()
+        return self._halfHeight
+
     def drawAt(self, z, axes, showLabels=False):  # pragma: no cover
         """ Draw element on plot with starting edge at 'z'.
 
@@ -972,22 +996,26 @@ class LensGraphic(MatrixGraphic):
                 if max(line._y) > maxRayHeight:
                     maxRayHeight = max(line._y)
 
-        halfHeight = self.displayHalfHeight(minSize=maxRayHeight)  # real units, i.e. data
+        self._halfHeight = self.displayHalfHeight(minSize=maxRayHeight)  # real units, i.e. data
 
         (xScaling, yScaling) = self.axesToDataScale(axes)
-        arrowHeadHeight = 2 * halfHeight * 0.08
+        arrowHeadHeight = 2 * self._halfHeight * 0.08
 
-        heightFactor = halfHeight * 2 / yScaling
+        heightFactor = self._halfHeight * 2 / yScaling
         arrowHeadWidth = xScaling * 0.008 * (heightFactor / 0.2) ** (3 / 4)
 
-        axes.arrow(z, 0, 0, halfHeight, width=arrowHeadWidth / 10, fc='k', ec='k',
+        axes.arrow(z, 0, 0, self._halfHeight, width=arrowHeadWidth / 10, fc='k', ec='k',
                    head_length=arrowHeadHeight, head_width=arrowHeadWidth, length_includes_head=True)
-        axes.arrow(z, 0, 0, -halfHeight, width=arrowHeadWidth / 10, fc='k', ec='k',
+        axes.arrow(z, 0, 0, -self._halfHeight, width=arrowHeadWidth / 10, fc='k', ec='k',
                    head_length=arrowHeadHeight, head_width=arrowHeadWidth, length_includes_head=True)
         self.drawCardinalPoints(z, axes)
 
 
 class SpaceGraphic(MatrixGraphic):
+    def __init__(self, matrix):
+        super(SpaceGraphic, self).__init__(matrix)
+        self._halfHeight = 0
+
     def drawAt(self, z, axes, showLabels=False):
         """This function draws nothing because free space is not visible. """
         return
@@ -1168,6 +1196,7 @@ class ApertureGraphic(MatrixGraphic):
 class MatrixGroupGraphic(MatrixGraphic):
     def __init__(self, matrixGroup: MatrixGroup):
         self.matrixGroup = matrixGroup
+        self._halfHeight = None
         super().__init__(matrixGroup)
 
     @property
@@ -1177,6 +1206,12 @@ class MatrixGroupGraphic(MatrixGraphic):
             L += element.L
         return L
 
+    @property
+    def halfHeight(self):
+        if self._halfHeight is None:
+            self._halfHeight = self.displayHalfHeight()
+        return self._halfHeight
+
     def drawAt(self, z, axes, showLabels=True):
         """ Draw each element of this group """
         for element in self.matrixGroup:
@@ -1220,13 +1255,41 @@ class MatrixGroupGraphic(MatrixGraphic):
                     labels[zStr] = label
             zElement += element.L
 
-        halfHeight = self.matrixGroup.largestDiameter() / 2
+        halfHeight = self.matrixGroup.largestDiameter / 2
         for zStr, label in labels.items():
             z = float(zStr)
             axes.annotate(label, xy=(z, 0.0), xytext=(z, -halfHeight * 0.5),
                           xycoords='data', fontsize=12,
                           ha='center', va='bottom')
 
+    def display(self):
+        fig, axes = plt.subplots(figsize=(10, 7))
+        self.drawAt(0, axes, showLabels=True)
+        self.drawAperture(0, axes)
+        self.drawPointsOfInterest(0, axes)
+        self.drawVertices(0, axes)
+        self.drawCardinalPoints(0, axes)
+        self.drawPrincipalPlanes(0, axes)
+        axes.set_ylim(-self.halfHeight * 1.6, self.halfHeight * 1.6)
+        self._showPlot()
+
+    def _showPlot(self):  # pragma: no cover
+        # internal, do not use
+        try:
+            plt.plot()
+            if sys.platform.startswith('win'):
+                plt.show()
+            else:
+                plt.draw()
+                while True:
+                    if plt.get_fignums():
+                        plt.pause(0.001)
+                    else:
+                        break
+
+        except KeyboardInterrupt:
+            plt.close()
+
 
 class AchromatDoubletLensGraphic(MatrixGroupGraphic):
     def drawAt(self, z, axes, showLabels=False):
@@ -1465,17 +1528,19 @@ class ObjectiveGraphic(MatrixGroupGraphic):
 
 class Graphic:
     def __new__(cls, element):
-        if type(element) is AchromatDoubletLens:
+        if type(element) is AchromatDoubletLens or issubclass(type(element), AchromatDoubletLens):
             return AchromatDoubletLensGraphic(element)
-        if type(element) is SingletLens:
+        if type(element) is SingletLens or issubclass(type(element), SingletLens):
             return SingletLensGraphic(element)
-        if issubclass(type(element), Objective):
+        if issubclass(type(element), Objective) or issubclass(type(element), Objective):
             return ObjectiveGraphic(element)
         if issubclass(type(element), MatrixGroup):
             return MatrixGroupGraphic(element)
 
         if type(element) is Lens:
             return LensGraphic(element)
+        if type(element) is ThickLens:
+            return ThickLensGraphic(element)
         if type(element) is Space:
             return SpaceGraphic(element)
         if type(element) is DielectricInterface:
diff --git a/raytracing/gaussianbeam.py b/raytracing/gaussianbeam.py
index f478cf2..9755bdb 100644
--- a/raytracing/gaussianbeam.py
+++ b/raytracing/gaussianbeam.py
@@ -144,10 +144,10 @@ class GaussianBeam(object):
             description += "w(z): {0:.3f}, ".format(self.w)
             description += "R(z): {0:.3f}, ".format(self.R)
             description += "z: {0:.3f}, ".format(self.z)
-            description += "λ: {0:.1f} nm\n".format(self.wavelength * 1e6)
+            description += "λ: {0:.1f} nm\n".format(self.wavelength)
             description += "zo: {0:.3f}, ".format(self.zo)
             description += "wo: {0:.3f}, ".format(self.wo)
             description += "wo position: {0:.3f} ".format(self.waistPosition)
             return description
         else:
-            return "Not valid complex radius of curvature"
+            return "Beam is not finite: q={0}".format(self.q)
diff --git a/raytracing/lasercavity.py b/raytracing/lasercavity.py
index 09ae7aa..6e085b4 100644
--- a/raytracing/lasercavity.py
+++ b/raytracing/lasercavity.py
@@ -48,19 +48,23 @@ class LaserCavity(LaserPath):
         round trip: you will need to duplicate elements in reverse
         and append them manually.
 
+        Knowing that q = (Aq + B)/(Cq + d), we get
+        Cq^2 + Dq = Aq + B, therefore:
+        Cq^2 + (D-A)q - B = 0
+
+        and q = - ((D-A) +- sqrt( (D-A)^2 - 4 C (-B)))/(2C)
+
         You will typically obtain two values, where only one is physical.
         There could be two modes in a laser with complex matrices (i.e.
         with gain), but this is not considered here.  See "Lasers" by Siegman.
         """
         if not self.hasPower:
             return None, None
-
         b = self.D - self.A
-        sqrtDelta = cmath.sqrt(b * b + 4.0 * self.B * self.C)
+        sqrtDelta = cmath.sqrt(b * b - 4.0 * self.C * (-self.B))
 
         q1 = (- b + sqrtDelta) / (2.0 * self.C)
         q2 = (- b - sqrtDelta) / (2.0 * self.C)
-
         return (GaussianBeam(q=q1), GaussianBeam(q=q2))
 
     def laserModes(self):
@@ -91,5 +95,8 @@ class LaserCavity(LaserPath):
             graph in the bottom half of the plot. (default=None)
 
         """
+        beams = self.laserModes()
+        if len(beams) == 0:
+            print("Cavity is not stable")
 
-        super(LaserCavity, self).display(inputBeams=self.laserModes())
+        super(LaserCavity, self).display(beams=beams)
diff --git a/raytracing/laserpath.py b/raytracing/laserpath.py
index fd2a04b..aa8ac42 100644
--- a/raytracing/laserpath.py
+++ b/raytracing/laserpath.py
@@ -53,24 +53,25 @@ class LaserPath(MatrixGroup):
         self.showPlanesAcrossPointsOfInterest = True
         super(LaserPath, self).__init__(elements=elements, label=label)
 
-    def display(self, inputBeams=None, comments=None):  # pragma: no cover
+    def display(self, beams=None, comments=None):  # pragma: no cover
         """ Display the optical system and trace the laser beam. 
         If comments are included they will be displayed on a
         graph in the bottom half of the plot.
 
         Parameters
         ----------
+        inputBeam : object of GaussianBeam class
         inputBeams : list of object of GaussianBeam class
             A list of Gaussian beams
         comments : string
             If comments are included they will be displayed on a graph in the bottom half of the plot. (default=None)
 
         """
-        if inputBeams is None:
-            inputBeams = [self.inputBeam]
+        if beams is None :
+            beams = [self.inputBeam]
 
         figure = Figure(opticalPath=self)
 
         figure.createFigure(title=self.label, comments=comments)
 
-        figure.displayGaussianBeam(inputBeams=inputBeams)
+        figure.displayGaussianBeam(beams=beams)
diff --git a/raytracing/matrix.py b/raytracing/matrix.py
index 2f8ff11..4c2a801 100644
--- a/raytracing/matrix.py
+++ b/raytracing/matrix.py
@@ -1364,6 +1364,10 @@ class Matrix(object):
             halfHeight = self.apertureDiameter / 2.0  # real half height
         return halfHeight
 
+    def display(self):
+        from .figure import Graphic
+        return Graphic(self).display()
+
     def __str__(self):
         """ String description that allows the use of print(Matrix())
 
diff --git a/raytracing/specialtylenses.py b/raytracing/specialtylenses.py
index 102815e..43dd49d 100644
--- a/raytracing/specialtylenses.py
+++ b/raytracing/specialtylenses.py
@@ -137,7 +137,6 @@ class AchromatDoubletLens(MatrixGroup):
         (f1, f2) = self.focusPositions(z)
         return [{'z': f1, 'label': '$F_f$'}, {'z': f2, 'label': '$F_b$'}]
 
-
 class SingletLens(MatrixGroup):
     """
     General singlet lens with an effective focal length of f, back focal


Test patch:
diff --git a/raytracing/tests/testsFigure.py b/raytracing/tests/testsFigure.py
index 0c29ffd..2ea8f0a 100644
--- a/raytracing/tests/testsFigure.py
+++ b/raytracing/tests/testsFigure.py
@@ -5,29 +5,34 @@ from raytracing import *
 
 
 class TestFigure(unittest.TestCase):
-    def testDisplayRangeWithFiniteLens(self):
+    @patch('raytracing.Figure._showPlot')
+    def testDisplayRangeWithFiniteLens(self, mock):
         path = ImagingPath()  # default objectHeight is 10
         path.append(Space(d=10))
         path.append(Lens(f=5, diameter=20))
+        path.display()
 
         largestDiameter = 20
 
         self.assertEqual(path.figure.displayRange(), largestDiameter)
 
-    def testDisplayRangeImageOutOfView(self):
+    @patch('raytracing.Figure._showPlot')
+    def testDisplayRangeImageOutOfView(self, mock):
         path = ImagingPath()
         path.append(Space(2))
         path.append(CurvedMirror(-5, 10))
+        path.display()
 
-        self.assertAlmostEqual(path.figure.displayRange(), 20)
+        self.assertAlmostEqual(path.figure.displayRange(), 10)
 
         path.objectHeight = 1
+        path.display()
         self.assertEqual(path.figure.displayRange(), 10)
 
     def testDisplayRangeWithEmptyPath(self):
         path = ImagingPath()
 
-        largestDiameter = path.objectHeight * 2
+        largestDiameter = path.objectHeight
 
         self.assertEqual(path.figure.displayRange(), largestDiameter)
 
@@ -70,19 +75,8 @@ class TestFigureAxesToDataScale(unittest.TestCase):
 
         xScaling, yScaling = figure.axesToDataScale()
 
-        self.assertEqual(xScaling, 1)
-        self.assertEqual(yScaling, 1)
-
-    def testWithForcedScale(self):
-        figure = Figure(ImagingPath())
-        figure.createFigure()
-        figure.axes.set_xlim(-10, 10)
-        figure.axes.set_ylim(-5, 5)
-
-        xScaling, yScaling = figure.axesToDataScale()
-
-        self.assertEqual(xScaling, 20)
-        self.assertEqual(yScaling, 10)
+        self.assertEqual(xScaling, 0)
+        self.assertEqual(yScaling, 10 * 1.6)
 
     @unittest.skipIf(sys.platform == 'darwin',"FIXME: We hacked plt.show() on darwin to recover Ctrl-C")
     @patch('matplotlib.pyplot.show')
@@ -96,7 +90,7 @@ class TestFigureAxesToDataScale(unittest.TestCase):
 
         (xScaling, yScaling) = path.figure.axesToDataScale()
 
-        self.assertEqual(yScaling, path.figure.displayRange() * 1.5)
+        self.assertEqual(yScaling, path.figure.displayRange() * 1.6)
         self.assertEqual(xScaling, 20 * 1.1)
 
 
diff --git a/raytracing/tests/testsGaussian.py b/raytracing/tests/testsGaussian.py
index 66e7e5c..6d04013 100644
--- a/raytracing/tests/testsGaussian.py
+++ b/raytracing/tests/testsGaussian.py
@@ -114,7 +114,8 @@ class TestBeam(envtest.RaytracingTestCase):
 
     def testStrInvalidRadiusOfCurvature(self):
         beam = GaussianBeam(w=inf, R=1)
-        self.assertEqual(str(beam), "Not valid complex radius of curvature")
+        self.assertFalse(beam.isFinite)
+        self.assertEqual(str(beam), "Beam is not finite: q=(1+0j)")
 
 
 if __name__ == '__main__':