Repository: DCC-Lab/Polarization (language: python)

Issue:
Speed optimization
The code as it stands provide the right answers (similar to MATLAB code), but is abysmally slow. It does not matter, this can be optimized.

There are several things that can be implemented to improve the performance considerably.  As it stands, there are 3 lines of optimizations to follow:
1. **The most important penalty is the redundancy of several calculations**: because we made the choice to encapsulate the k-dependence in the `JonesMatrix` (i.e. the numerical values of certain matrices are not available until the `JonesVector` is known), then we often end up re-calculating a matrix for a given k over and over.  The obvious optimization here is to cache the `forward` and `backward` matrices for a given `k` and simply return the computed values if it is for the same `k`. This all boils down to implementing a caching mechanism in `computeMatrix` for `JonesMatrix`, but also for `MatrixProduct`.
2. The object penalty.  Going through objects for various calculations adds a performance penalty.  For instance, `@properties` are slower than variables, `__init__` takes time, etc... In some cases, this penalty is unacceptable. The best example is the *internal* calculation of the backscatterMatrix for a TissueLayer from all scatterers. It is better to explicitly compute a single matrix rather than take all the individual matrices from each scatterers.
3. Finally, there appear to be various performance enhancements possible with numpy versus other means of calculations but this looks less promizing.

This could be approached two ways:
1. Local caching mechanisms in objects
2. Global pre-flight optimization, pre computing everything.

Since the real problem is the redundancy of the calculations, I don't think the global strategy will blow the local strategy out of the water: the real gain is to compute everything just once, and the local/global strategies are two means to an end. I would much prefer a local optimization, but we may experiment with both.

Let's implement these changes individually in a unique branche for each optimization.
Let's always confirm the unittests are still passing, because the current branch gives results that are similar to the original MATLAB code.
PS Tissue Tests
Added quantitative tests for PS propagation and backscattering in tissue. 


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 5cd6377..d4334ea 100644
--- a/.gitignore
+++ b/.gitignore
@@ -124,3 +124,4 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+data/
\ No newline at end of file
diff --git a/polarization/jonesmatrix.py b/polarization/jonesmatrix.py
index e768b0b..c98a153 100644
--- a/polarization/jonesmatrix.py
+++ b/polarization/jonesmatrix.py
@@ -11,13 +11,14 @@ class JonesMatrix:
     It is a complex 2 x 2 matrix that can transform a Jones vector
     (representing Ex and Ey).
     """
+
     def __init__(self, A: complex = None,
-                       B: complex = None,
-                       C: complex = None,
-                       D: complex = None, 
-                       m=None, 
-                       physicalLength: float = 0,
-                       orientation: float = 0):
+                 B: complex = None,
+                 C: complex = None,
+                 D: complex = None,
+                 m=None,
+                 physicalLength: float = 0,
+                 orientation: float = 0):
         """ We may initialize the matrix with all four elements or with a numpy
         array. It is mandatory to provide a length for the element, and it is
         assumed that the basis for the matrix is x̂ and ŷ.
@@ -27,22 +28,60 @@ def __init__(self, A: complex = None,
         to avoid problems. For now, it is x and y.
         """
         if m is not None:
-            self.mOriginal = array(m)
+            self.mOriginal = m
         elif A is not None and D is not None and C is not None and D is not None:
-            self.mOriginal = array([[A,B],[C,D]])
+            self.mOriginal = [A, B, C, D]
         else:
             # Obviously, the subclass will compute when the time comes
             # See BirefringenceMaterial() for an example
             self.mOriginal = None
 
-        self.orientation = orientation
         self.L = physicalLength
 
+        self.orientation = orientation
+
         """ The basis vector for x and y. For now this is not really
         modifiable.  b1 x b2 = b3, direction of propagation """
-        self.b1 = Vector(1,0,0) # x̂
-        self.b2 = Vector(0,1,0) # ŷ
-        self.b3 = Vector(0,0,1) # ẑ
+        self.b1 = xHat  # Vector(1,0,0) # x̂
+        self.b2 = yHat  # Vector(0,1,0) # ŷ
+        self.b3 = zHat  # Vector(0,0,1) # ẑ
+
+    @property
+    def mOriented(self):
+        value = self.orientation
+        c = complex(cos(value))
+        s = complex(sin(value))
+        c2 = c * c
+        s2 = s * s
+        cs = c * s
+
+        A = self.mOriginal[0]
+        B = self.mOriginal[1]
+        C = self.mOriginal[2]
+        D = self.mOriginal[3]
+
+        if B == 0 and C == 0:
+            Ap = c2 * A + D * s2
+            Bp = A * cs - cs * D
+            Cp = Bp
+            Dp = c2 * D + A * s2
+        else:
+            Ap = c2 * A - C * cs - cs * B + D * s2
+            Bp = c2 * B + A * cs - cs * D - C * s2
+            Cp = c2 * C - D * cs + cs * A - B * s2
+            Dp = c2 * D + B * cs + cs * C + A * s2
+
+        return [Ap, Bp, Cp, Dp]
+
+    def computePythonMatrix(self, k=None, l=None, backward=bool):
+        if self.mOriginal is None:
+            # This is the signal that the matrix depends on k.
+            # subclasses will calculate m as a function of k when needed
+            # by overriding this function.
+            raise ValueError('This matrix {0} appears to be wavelength-dependent. \
+You cannot obtain the values without providing a wavevector k or the matrix itself.'.format(type(self)))
+
+        return self.mOriented
 
     def computeMatrix(self, k=None, l=None, backward=bool):
         if self.mOriginal is None:
@@ -52,12 +91,9 @@ def computeMatrix(self, k=None, l=None, backward=bool):
             raise ValueError('This matrix {0} appears to be wavelength-dependent. \
 You cannot obtain the values without providing a wavevector k or the matrix itself.'.format(type(self)))
 
-        theta = self.orientation
-        rotMatrix = array([[cos(theta),sin(theta)],[-sin(theta), cos(theta)]],dtype='complex')
-        invRotMatrix = array([[cos(theta),-sin(theta)],[sin(theta), cos(theta)]],dtype='complex')
+        A, B, C, D = self.computePythonMatrix()
+        return array([[A, B], [C, D]])
 
-        return matmul(invRotMatrix, matmul(self.mOriginal, rotMatrix))
-    
     def backward(self):
         """ Return a matrix for a JonesVector propagating in the opposite
         direction. Flip the direction of propagation that is assumed when  getting the
@@ -72,50 +108,49 @@ def backward(self):
         or Faraday rotator.
         FIXME: I am not sure about the orientation.
         """
-
-        backward = JonesMatrix(m=self.m.T,
+        A, B, C, D = self.m
+        backward = JonesMatrix(m=[A, C, B, D],
                                physicalLength=self.L,
-                               orientation=0)
+                               orientation=0)  # orientation already included inside self.m
         backward.b3 = -backward.b3
         backward.b2 = -backward.b2
-        backward.orientation = -backward.orientation
         return backward
 
     @property
     def m(self):
-        return self.computeMatrix()
-    
+        return self.computePythonMatrix()
+
     @property
     def A(self):
-        return self.m[0,0]
+        return self.m[0]
 
     @property
     def B(self):
-        return self.m[0,1]
+        return self.m[1]
 
     @property
     def C(self):
-        return self.m[1,0]
+        return self.m[2]
 
     @property
     def D(self):
-        return self.m[1,1]
+        return self.m[3]
 
     @property
     def determinant(self):
-        return det(self.m)
-    
+        return det(self.computeMatrix())
+
     @property
     def isBirefringent(self) -> bool:
         """ Returns True if it is birefringent.  See birefringence."""
 
         phi1, phi2, e1, e2 = self.birefringence
-        if isNotZero(phi1-phi2, epsilon=1e-7):
+        if isNotZero(phi1 - phi2, epsilon=1e-7):
             return True
         return False
 
     @property
-    def birefringence(self) :
+    def birefringence(self):
         """ Returns the birefringence parameters of this element. To determine
         these parameters, one must have a diagonalized version of this matrix
         and then determine if the diagonal element have a phase difference.
@@ -156,7 +191,7 @@ def eigens(self):
         the eigenvectors real by transferring a j factor to the eigenvalue if
         the eigenvector is imaginary """
 
-        w, v = eig(self.m)
+        w, v = eig(self.computeMatrix())
 
         e1 = realIfPossible(v[0])
         e2 = realIfPossible(v[1])
@@ -179,8 +214,8 @@ def diattenuation(self) -> complex:
     @property
     def retardance(self) -> float:
         phi1, phi2, e1, e2 = self.birefringence
-        return phi1-phi2
-    
+        return phi1 - phi2
+
     def __mul__(self, rightSide):
         """Operator overloading allowing easy-to-read matrix multiplication
         with other `JonesMatrix`, with a `JonesVector`.
@@ -190,14 +225,14 @@ def __mul__(self, rightSide):
         with rOut = M1*r
 
         """
-        if isinstance(rightSide, JonesMatrix):
+        if isinstance(rightSide, JonesVector):
+            return self.mul_vector(rightSide)
+        elif isinstance(rightSide, JonesMatrix):
             return self.mul_matrix(rightSide)
         elif isinstance(rightSide, MatrixProduct):
             product = MatrixProduct(matrices=rightSide.matrices)
             product.append(self)
             return product
-        elif isinstance(rightSide, JonesVector):
-            return self.mul_vector(rightSide)
         elif isinstance(rightSide, number_types):
             return self.mul_number(rightSide)
         else:
@@ -213,6 +248,7 @@ def __rmul__(self, leftSide):
             raise TypeError(
                 "Unrecognized left side element in multiply: '{0}'\
                  cannot be multiplied by a JonesMatrix".format(leftSide))
+
     def mul_matrix(self, rightSideMatrix: 'JonesMatrix'):
         r""" This function is used to combine two elements into a single matrix.
 
@@ -239,14 +275,14 @@ def mul_matrix(self, rightSideMatrix: 'JonesMatrix'):
 
         try:
             theMatrix = self.computeMatrix()
-            product = JonesMatrix(m=matmul(theMatrix, rightSideMatrix.m), physicalLength=self.L + rightSideMatrix.L)
+            m = matmul(theMatrix, rightSideMatrix.computeMatrix())
+            product = JonesMatrix(m=m.flatten(), physicalLength=self.L + rightSideMatrix.L)
             return product
         except ValueError as err:
             # There is no possible numerical value at this point. Let's return an
             # object the holds the elements together to be multiplied later
             # when we know the JonesVector
-            return MatrixProduct( [rightSideMatrix, self])
-
+            return MatrixProduct([rightSideMatrix, self])
 
     def mul_vector(self, rightSideVector):
         r"""This function does the multiplication of a vector by a matrix.
@@ -264,33 +300,30 @@ def mul_vector(self, rightSideVector):
 
         """
 
-        outputVector = JonesVector()
         # We obtain the matrix specific to this JonesVector
-        m = self.computeMatrix(k = rightSideVector.k)
-        
+        m = self.computePythonMatrix(k=rightSideVector.k)
+
         # Is the matrix in the appropriate direction?
         # For now, print a warning.
-        if rightSideVector.b3 != self.b3:
-            print("Warning: the matrix is set explicitly set up for propagation along the opposite direction of the JonesVector")
+        # if rightSideVector.b3 != self.b3:
+        #     print("Warning: the matrix is set explicitly set up for propagation along the opposite direction of the JonesVector")
 
-        direction = zHat.dot(rightSideVector.b3) # +1 or -1
+        direction = rightSideVector.b3.z  # zHat.dot(rightSideVector.b3) # +1 or -1
 
-        outputVector.Ex = m[0,0] * rightSideVector.Ex + m[0,1] * rightSideVector.Ey
-        outputVector.Ey = m[1,0] * rightSideVector.Ex + m[1,1] * rightSideVector.Ey
-        outputVector.z = rightSideVector.z + direction*self.L
-        outputVector.k = rightSideVector.k
-        return outputVector
+        E1 = m[0] * rightSideVector.E1 + m[1] * rightSideVector.E2
+        E2 = m[2] * rightSideVector.E1 + m[3] * rightSideVector.E2
+        return JonesVector(Ex=E1, Ey=E2, z=rightSideVector.z + direction * self.L, k=rightSideVector.k)
 
     def mul_number(self, n):
         """ Multiply a Jones matrix by a number."""
-        return JonesMatrix(self.A*n, self.B*n, self.C*n, self.D*n, physicalLength=self.L)
+        return JonesMatrix(self.A * n, self.B * n, self.C * n, self.D * n, physicalLength=self.L)
 
     def rotatedBy(self, theta):
         """ We return a rotated copy of the optical element of the matrix by theta. 
         For instance, a theta rotation of a horizontal polarizer will be a polarizer
         aligned at theta. """
 
-        return Rotation(theta = -theta)*self*Rotation(theta = theta)
+        return Rotation(theta=-theta) * self * Rotation(theta=theta)
 
     def rotateBasisBy(self, theta):
         """ We rotate the basis of the matrix by theta. This is different 
@@ -309,7 +342,7 @@ def setValue(self, name, value):
     def value(self, name):
         return getattr(self, name)
 
-    def show(self, input, xObj, xProperty, xRange, yObj, yProperty): # pragma: no cover
+    def show(self, input, xObj, xProperty, xRange, yObj, yProperty):  # pragma: no cover
         if xObj is None:
             xObj = self
 
@@ -322,51 +355,51 @@ def show(self, input, xObj, xProperty, xRange, yObj, yProperty): # pragma: no co
             else:
                 xObj.setValue(xProperty, value)
 
-            vOut = self*input
+            vOut = self * input
 
-            x.append( xObj.value(xProperty))
+            x.append(xObj.value(xProperty))
             if yObj is None:
-                y.append( vOut.value(yProperty))
+                y.append(vOut.value(yProperty))
             else:
-                y.append( yObj.value(yProperty))
-        
+                y.append(yObj.value(yProperty))
+
         plt.title("{0} versus {1}".format(yProperty, yProperty))
         plt.xlabel("{0}".format(xProperty))
         plt.ylabel("{0}".format(yProperty))
-        plt.plot(x,y,'ko')
+        plt.plot(x, y, 'ko')
         plt.show()
 
-    def showOrientationPlot(self, input:JonesVector): # pragma: no cover
+    def showOrientationPlot(self, input: JonesVector):  # pragma: no cover
         x = []
         y = []
-        for theta in range(0,190,10):
+        for theta in range(0, 190, 10):
             theMatrix = JonesMatrix(m=self.m)
-            theMatrix.orientation = theta*radPerDeg
+            theMatrix.orientation = theta * radPerDeg
 
-            vOut = theMatrix*input
+            vOut = theMatrix * input
             x.append(theta)
             y.append(vOut.intensity)
 
         plt.title("Intensity versus orientation of element")
         plt.xlabel(r"Rotation $\theta$ of element [°]")
         plt.ylabel(r"Intensity [arb. unit]")
-        plt.plot(x,y,'ko')
+        plt.plot(x, y, 'ko')
         plt.show()
 
-    def showPolarizationPlot(self): # pragma: no cover
+    def showPolarizationPlot(self):  # pragma: no cover
         x = []
         y = []
-        for theta in range(0,190,10):
-            vIn = JonesVector(Ex=cos(theta*radPerDeg), Ey=sin(theta*radPerDeg))
-            vOut = self*vIn
+        for theta in range(0, 190, 10):
+            vIn = JonesVector(Ex=cos(theta * radPerDeg), Ey=sin(theta * radPerDeg))
+            vOut = self * vIn
             x.append(theta)
             y.append(vOut.intensity)
 
         plt.title("Intensity versus input polarization orientation (0° is horizontal)")
         plt.xlabel(r"Input polarization orientation [°]")
-        plt.xlim(0,180)
+        plt.xlim(0, 180)
         plt.ylabel(r"Intensity [arb. unit]")
-        plt.plot(x,y,'ko')
+        plt.plot(x, y, 'ko')
         plt.show()
 
     def __str__(self):
@@ -375,38 +408,46 @@ def __str__(self):
 
 class HorizontalPolarizer(JonesMatrix):
     def __init__(self):
-        JonesMatrix.__init__(self, A=1,B=0,C=0,D=0,physicalLength=0)        
+        JonesMatrix.__init__(self, A=1, B=0, C=0, D=0, physicalLength=0)
+
 
 class LinearPolarizer(JonesMatrix):
     def __init__(self, theta):
         c = cos(theta)
         s = sin(theta)
-        JonesMatrix.__init__(self, A=c*c,B=c*s,C=c*s,D=s*s,physicalLength=0)
+        JonesMatrix.__init__(self, A=c * c, B=c * s, C=c * s, D=s * s, physicalLength=0)
+
 
 class VerticalPolarizer(JonesMatrix):
     def __init__(self):
-        JonesMatrix.__init__(self, A=0,B=0,C=0,D=1,physicalLength=0)        
+        JonesMatrix.__init__(self, A=0, B=0, C=0, D=1, physicalLength=0)
+
 
 class Plus45Polarizer(JonesMatrix):
     def __init__(self):
-        JonesMatrix.__init__(self, A=0.5,B=0.5,C=0.5,D=0.5,physicalLength=0)        
+        JonesMatrix.__init__(self, A=0.5, B=0.5, C=0.5, D=0.5, physicalLength=0)
+
 
 class Minus45Polarizer(JonesMatrix):
     def __init__(self):
-        JonesMatrix.__init__(self, A=0.5,B=-0.5,C=-0.5,D=0.5,physicalLength=0)        
+        JonesMatrix.__init__(self, A=0.5, B=-0.5, C=-0.5, D=0.5, physicalLength=0)
+
 
 class RightCircularPolarizer(JonesMatrix):
     def __init__(self):
-        JonesMatrix.__init__(self, A=0.5,B=0.5j,C=-0.5j,D=0.5,physicalLength=0)        
+        JonesMatrix.__init__(self, A=0.5, B=0.5j, C=-0.5j, D=0.5, physicalLength=0)
+
 
 class LeftCircularPolarizer(JonesMatrix):
     def __init__(self):
-        JonesMatrix.__init__(self, A=0.5,B=-0.5j,C=0.5j,D=0.5,physicalLength=0)        
+        JonesMatrix.__init__(self, A=0.5, B=-0.5j, C=0.5j, D=0.5, physicalLength=0)
+
 
 class Rotation(JonesMatrix):
     def __init__(self, theta):
         JonesMatrix.__init__(self, A=cos(theta), B=sin(theta), C=-sin(theta), D=cos(theta), physicalLength=0)
 
+
 class PhaseRetarder(JonesMatrix):
     def __init__(self, delta=None, phiX=None, phiY=None, physicalLength=0):
         if delta is not None:
@@ -414,25 +455,33 @@ def __init__(self, delta=None, phiX=None, phiY=None, physicalLength=0):
         else:
             JonesMatrix.__init__(self, A=exp(1j * phiX), B=0, C=0, D=exp(1j * phiY), physicalLength=physicalLength)
 
+
 class BirefringentMaterial(JonesMatrix):
-    def __init__(self, deltaIndex:float, fastAxisOrientation, physicalLength=0):
+    def __init__(self, deltaIndex: float, fastAxisOrientation, physicalLength=0):
         """ The fast axis is the X axis when fastAxisOrientation = 0"""
-        JonesMatrix.__init__(self, A=None, B=None, C=None, D=None, physicalLength=physicalLength, orientation=fastAxisOrientation)
+        JonesMatrix.__init__(self, A=None, B=None, C=None, D=None, physicalLength=physicalLength,
+                             orientation=fastAxisOrientation)
         self.deltaIndex = deltaIndex
         self.isBackward = False
+        c = complex(cos(self.orientation))
+        s = complex(sin(self.orientation))
+        self.c2 = c * c
+        self.s2 = s * s
+        self.cs = c * s
 
     def computeMatrix(self, k=None):
+        Ap, Bp, Bp, Dp = self.computePythonMatrix(k)
+        return array([[Ap, Bp], [Bp, Dp]])
+
+    def computePythonMatrix(self, k=None):
         if k is not None:
-            phi = k * self.L
-            deltaPhi = k * self.deltaIndex * self.L
-            explicit = JonesMatrix(A=exp(1j * phi), B=0, C=0, D=exp(1j * (deltaPhi + phi)), physicalLength=self.L)
-            explicit.orientation = self.orientation
-            if self.isBackward:
-                explicit = JonesMatrix(m=explicit.m.T, physicalLength=self.L,
-                                       orientation=0)
-                explicit.b3 = -explicit.b3
-                explicit.b2 = -explicit.b2
-            return explicit.computeMatrix()
+            jPhi = 1j * self.L * k
+            A = exp(jPhi)
+            D = exp(jPhi * (1 + self.deltaIndex))
+            Ap = self.c2 * A + D * self.s2
+            Bp = A * self.cs - self.cs * D
+            Dp = self.c2 * D + A * self.s2
+            return [Ap, Bp, Bp, Dp]
         else:
             raise ValueError("You must provide k for this matrix")
 
@@ -440,9 +489,11 @@ def backward(self):
         backward = BirefringentMaterial(deltaIndex=self.deltaIndex,
                                         fastAxisOrientation=self.orientation,
                                         physicalLength=self.L)
-        backward.isBackward = True
+        backward.b3 = -self.b3
+        backward.b2 = -self.b2
         return backward
 
+
 class Vacuum(JonesMatrix):
     def __init__(self, physicalLength=0):
         """ The fast axis is the X axis when fastAxisOrientation = 0"""
@@ -450,38 +501,41 @@ def __init__(self, physicalLength=0):
         self.isBackward = False
 
     def computeMatrix(self, k=None):
+        Ap, Bp, Bp, Dp = self.computePythonMatrix(k)
+        return array([[Ap, Bp], [Bp, Dp]])
+
+    def computePythonMatrix(self, k=None):
         if k is not None:
-            explicit = JonesMatrix(A=exp(1j * k * self.L), B=0, C=0, D=exp(1j * k * self.L), physicalLength=self.L)
-            explicit.orientation = self.orientation
-            if self.isBackward:
-                explicit = JonesMatrix(m=explicit.m.T, physicalLength=self.L,
-                                       orientation=0)
-                explicit.b3 = -explicit.b3
-                explicit.b2 = -explicit.b2
-            return explicit.computeMatrix()
+            A = exp(1j * k * self.L)
+            return [A, 0, 0, A]
         else:
             raise ValueError("You must provide k for this matrix")
 
     def backward(self):
         backward = Vacuum(physicalLength=self.L)
-        backward.isBackward = True
+        backward.b3 = -self.b3
+        backward.b2 = -self.b2
         return backward
 
+
 class Diattenuator(JonesMatrix):
     def __init__(self, Tx, Ty, physicalLength=0):
         JonesMatrix.__init__(self, A=Tx, B=0, C=0, D=Ty, physicalLength=physicalLength)
-    
+
+
 class QWP(PhaseRetarder):
     def __init__(self, theta):
         # theta is fast axis with respect to x-axis
-        PhaseRetarder.__init__(self, delta=-pi / 2) # Ex is advanced by pi/2, x is fast
+        PhaseRetarder.__init__(self, delta=-pi / 2)  # Ex is advanced by pi/2, x is fast
         self.orientation = theta
 
+
 class HWP(PhaseRetarder):
     def __init__(self, theta):
-        PhaseRetarder.__init__(self, delta=-pi) # Ex is advanced by pi, x is fast
+        PhaseRetarder.__init__(self, delta=-pi)  # Ex is advanced by pi, x is fast
         self.orientation = theta
 
+
 class PockelsCell(JonesMatrix):
     def __init__(self, halfwaveVoltage, length):
         self.halfwaveVoltage = halfwaveVoltage
@@ -490,32 +544,32 @@ def __init__(self, halfwaveVoltage, length):
 
     @property
     def m(self):
-        cell = PhaseRetarder(delta=self.voltage/self.halfwaveVoltage*pi)
+        cell = PhaseRetarder(delta=self.voltage / self.halfwaveVoltage * pi)
         cell.orientation = self.orientation
         return cell.m
 
-    def showVoltagePlot(self): # pragma: no cover
+    def showVoltagePlot(self):  # pragma: no cover
         fig, axs = plt.subplots(1, sharex=True)
-        fig.suptitle("Pockels cell at {0:.1f}° with horizontal incident polarization".format(self.orientation*degPerRad))
+        fig.suptitle(
+            "Pockels cell at {0:.1f}° with horizontal incident polarization".format(self.orientation * degPerRad))
 
-        voltages = list(range(0,self.halfwaveVoltage+10,10))
+        voltages = list(range(0, self.halfwaveVoltage + 10, 10))
         yParallel = []
         yCrossed = []
         for voltage in voltages:
             self.voltage = voltage
 
             vIn = JonesVector.horizontal()
-            vOut = HorizontalPolarizer()*self*vIn
+            vOut = HorizontalPolarizer() * self * vIn
             yParallel.append(vOut.intensity)
 
-            vOut = VerticalPolarizer()*self*vIn
+            vOut = VerticalPolarizer() * self * vIn
             yCrossed.append(vOut.intensity)
 
-
         # axs[0].xlim(0, self.halfwaveVoltage)
         # fig.ylabel(r"Intensity [arb. unit]")
-        axs.plot(voltages,yParallel,'k|',label="Between parallel polarizers")
-        axs.plot(voltages,yCrossed,'k+',label="Between crossed polarizers")
+        axs.plot(voltages, yParallel, 'k|', label="Between parallel polarizers")
+        axs.plot(voltages, yCrossed, 'k+', label="Between crossed polarizers")
         axs.legend()
         plt.show()
 
@@ -577,9 +631,8 @@ def mul_vector(self, vector):
         and multiply the rightmost matrix by the JonesVector , and if that  matrix
         requires the vector k, it will request it in mul_vector in order to calculate
         the numerical value of the matrix. """
-        outputVector = JonesVector(Ex=vector.Ex, Ey=vector.Ey, k=vector.k, z=vector.z)
+        outputVector = JonesVector(Ex=vector.E1, Ey=vector.E2, k=vector.k, z=vector.z)
         for m in self.matrices:
-            outputVector = m*outputVector
+            outputVector.transformBy(m)
 
         return outputVector
-
diff --git a/polarization/jonesvector.py b/polarization/jonesvector.py
index 7c6c868..4577a99 100644
--- a/polarization/jonesvector.py
+++ b/polarization/jonesvector.py
@@ -22,9 +22,9 @@ def __init__(self, Ex: complex = 0.0, Ey: complex = 0.0, wavelength:float = None
         if k is not None and wavelength is not None:
             raise ValueError('Provide either one of wavelength or k, but not both')
         elif k is not None:
-            self.k = k
+            self.k = complex(k)
         elif wavelength is not None:
-            self.k = 2*np.pi/wavelength
+            self.k = complex(2*np.pi/wavelength)
         else:
             # If at any point in calculation k is needed, it will fail
             # This is the expected behaviour: k must be set explicitly
@@ -33,9 +33,9 @@ def __init__(self, Ex: complex = 0.0, Ey: complex = 0.0, wavelength:float = None
 
         """ The basis vector for E1 and E2. We settled for b1, b2 and b3. For
         now this is not modifiable. """
-        self.b1 = Vector(1,0,0) # x̂ for E1
-        self.b2 = Vector(0,1,0) # ŷ for E2
-        self.b3 = Vector(0,0,1) # ẑ direction propagation b1 x b2 == b3
+        self.b1 = xHat # x̂ for E1
+        self.b2 = yHat # ŷ for E2
+        self.b3 = zHat # ẑ direction propagation b1 x b2 == b3
 
     @property
     def Ex(self):
@@ -106,7 +106,13 @@ def __add__(self, rhs):
         # if self.z != rhs.z:
         #     print("Warning: addition of two Jonesvectors from two different z: {0} and {1}".format(self.z, rhs.z))
 
-        return JonesVector(Ex=self.Ex+rhs.Ex, Ey=self.Ey+rhs.Ey, k=self.k, z=self.z)
+        return JonesVector(Ex=self.E1+rhs.E1, Ey=self.E2+rhs.E2, k=self.k, z=self.z)
+
+    def transformBy(self, rightSide):
+        m = rightSide.computePythonMatrix(self.k)
+        E1 = self.E1
+        self.E1 = E1*m[0] + self.E2*m[1]
+        self.E2 = E1*m[2] + self.E2*m[3]
 
     def __mul__(self, rightSide):
         if isinstance(rightSide, number_types):
@@ -129,7 +135,7 @@ def mul_number(self, n):
         return JonesVector(Ex=self.Ex * n, Ey=self.Ey * n, k=self.k, z=self.z)
 
     def copy(self):
-        return JonesVector(Ex=self.Ex, Ey=self.Ey, k=self.k, z=self.z)
+        return JonesVector(Ex=self.E1, Ey=self.E2, k=self.k, z=self.z)
 
     def reflect(self):
         
diff --git a/polarization/pulse.py b/polarization/pulse.py
index 3628d4e..da6fa14 100644
--- a/polarization/pulse.py
+++ b/polarization/pulse.py
@@ -142,7 +142,16 @@ def __init__(self, pulses):
         self.pulses = pulses
 
     @property
-    def fringes(self):  # or tomogram?
+    def k(self):
+        k = self.pulses[0].k
+        for pulse in self.pulses[1:]:
+            if pulse.k != k:
+                raise Exception("Undefined 'k' for PulseCollection since its pulses don't have "
+                                "the same 'k' distribution. ")
+        return k
+
+    @property
+    def fringes(self):
         """ Fringes of shape (2xN_states, width, resolution) """
         out = []
         for pulse in self.pulses:
@@ -169,12 +178,17 @@ def isExpanded(self):
                 return True
         return False
 
-    def display(self):
-        fig, axes = plt.subplots(1, 4)
+    def display(self, tag=""):
+        fig, axes = plt.subplots(1, 4, figsize=(14, 10))
         for i, ax in enumerate(axes):
             ax.imshow(np.transpose(self.intensity[i]), aspect='auto', cmap='gray')
+        plt.tight_layout()
+        # plt.savefig("test{}.png".format(tag))
         plt.show()
 
+    def save(self, filePath):
+        np.save(filePath, self.fringes)
+
     @classmethod
     def dualInputStates(cls, centerWavelength, wavelengthBandwidth, resolution=512):
         p1 = Pulse.horizontal(centerWavelength, wavelengthBandwidth, resolution=resolution)
diff --git a/polarization/tissue.py b/polarization/tissue.py
index dde0123..1877438 100644
--- a/polarization/tissue.py
+++ b/polarization/tissue.py
@@ -1,6 +1,5 @@
 from .tissueStack import *
 from .pulse import *
-from copy import deepcopy
 from typing import Union
 import numpy as np
 
@@ -17,6 +16,11 @@ def __init__(self, stacks: List[TissueStack], height=3000, width=200, depth=1):
         self.width = width
         self.depth = depth
 
+        self._scattDensity = None
+        self._opticAxis = None
+        self._apparentOpticAxis = None
+        self._birefringence = None
+
     def scan(self, pulse: Union[Pulse, PulseCollection], verbose=False):
         if verbose:
             def v_print(*args, **kwargs):
@@ -42,11 +46,14 @@ def _scanPulse(self, pulse, v_print) -> PulseArray:
     def _scanPulseCollection(self, pulses, v_print) -> PulseCollection:
         assert not pulses.isExpanded, "Cannot scan a PulseCollection that was already scanned. "
 
-        pulsesOut = []
-        for i, pulse in enumerate(pulses):
-            v_print("Pulse {}/{}".format(i+1, len(pulses)))
-            pulsesOut.append(self._scanPulse(pulse, v_print))
-        return PulseCollection(pulsesOut)
+        pulsesBScan = [[] for _ in pulses]
+        for a in range(self.width):
+            v_print(" .Stack {}/{}".format(a+1, self.width))
+            pulsesALine = self.stacks[a].backscatterMany(pulses)
+            for p in range(len(pulses)):
+                pulsesBScan[p].append(pulsesALine[p])
+
+        return PulseCollection(pulses=[PulseArray(bScan) for bScan in pulsesBScan])
 
     def __iter__(self):
         return iter(self.stacks)
@@ -54,9 +61,71 @@ def __iter__(self):
     def __len__(self):
         return len(self.stacks)
 
-    def display(self):
+    @property
+    def scattDensity(self):
+        if self._scattDensity is None:
+            self._scattDensity = self.stacks[0].scattDensity[:, None]
+            for stack in self.stacks[1:]:
+                self._scattDensity = np.concatenate([self._scattDensity, stack.scattDensity[:, None]], axis=1)
+        return self._scattDensity
+
+    @property
+    def opticAxis(self):
+        if self._opticAxis is None:
+            self._opticAxis = self.stacks[0].opticAxis[:, :, None]
+            for stack in self.stacks[1:]:
+                self._opticAxis = np.concatenate([self._opticAxis, stack.opticAxis[:, :, None]], axis=2)
+        return self._opticAxis
+
+    @property
+    def apparentOpticAxis(self):
+        if self._apparentOpticAxis is None:
+            self._apparentOpticAxis = self.stacks[0].apparentOpticAxis[:, :, None]
+            for stack in self.stacks[1:]:
+                self._apparentOpticAxis = np.concatenate([self._apparentOpticAxis, stack.apparentOpticAxis[:, :, None]], axis=2)
+        return self._apparentOpticAxis
+
+    @property
+    def birefringence(self):
+        if self._birefringence is None:
+            self._birefringence = self.stacks[0].birefringence[:, None]
+            for stack in self.stacks[1:]:
+                self._birefringence = np.concatenate([self._birefringence, stack.birefringence[:, None]], axis=1)
+        return self._birefringence
+
+    def display(self, title=None):
         """ Display all layer stacks and their properties. """
-        pass
+        fig, axes = plt.subplots(1, 5, figsize=(15, 5), sharey="all", sharex="all")
+        self._displayToAxes(axes)
+
+        fig.suptitle(title, fontsize=14)
+        plt.show()
+
+    def _displayToAxes(self, axes):
+        axes[0].imshow(self.scattDensity, aspect='auto', vmin=0, vmax=20, interpolation='none')
+        axes[0].set_title("Scattering Density")
+        axes[1].imshow(self.opticAxis[0], aspect='auto', vmin=-1, vmax=1, cmap='bwr', interpolation='none')
+        axes[1].set_title("OA1")
+        axes[2].imshow(self.opticAxis[1], aspect='auto', vmin=-1, vmax=1, cmap='bwr', interpolation='none')
+        axes[2].set_title("OA2")
+        axes[3].imshow(self.opticAxis[2], aspect='auto', vmin=-1, vmax=1, cmap='bwr', interpolation='none')
+        axes[3].set_title("OA3")
+        axes[4].imshow(self.birefringence, aspect='auto', vmin=0, vmax=0.0042, interpolation='none')
+        # todo: link vmax to sim layer.max_dn...
+        axes[4].set_title("Birefringence")
+
+    def save(self, filePath, reSample: int = None):
+
+        sampleData = np.concatenate([self.scattDensity[None, :, :],
+                                     # self.apparentOpticAxis,  # todo
+                                     self.opticAxis,
+                                     self.birefringence[None, :, :]])
+
+        if reSample is not None:
+            picks = np.linspace(0, sampleData.shape[1]-1, num=reSample, dtype=np.int)
+            sampleData = np.take(sampleData, picks, axis=1)
+
+        np.save(filePath, np.moveaxis(sampleData, 2, 1))
 
 
 class RandomTissue2D(Tissue):
@@ -75,7 +144,7 @@ def __init__(self, height=3000, width=200,
 
         self.flat = flat
         self._layerSizeMap = None
-        self.referenceStack = referenceStack
+        self.referenceStack: TissueStack = referenceStack
 
         self.generateMap()
         self.generateStacks()
@@ -112,12 +181,16 @@ def referenceStack(self, stack):
 
     def _stackOf(self, layerSizes):
         layers = []
-        for thickness, layer in zip(layerSizes[1:], deepcopy(self.referenceStack.layers)):
-            layer.thickness = thickness
-            layer.resetScatterers()
-            layers.append(layer)
-
-        return TissueStack(offset=layerSizes[0], layers=layers)
+        currentHeight = layerSizes[0]
+        for thickness, layer in zip(layerSizes[1:], self.referenceStack.layers):
+            if currentHeight >= self.height:
+                continue
+            elif currentHeight + thickness > self.height:
+                thickness = self.height - currentHeight
+            layers.append(layer.copy(thickness=thickness))
+            currentHeight += thickness
+        stack = TissueStack(offset=layerSizes[0], layers=layers, height=self.height)
+        return stack
 
     def generateStacks(self):
         for w in range(self.width):
diff --git a/polarization/tissueLayer.py b/polarization/tissueLayer.py
index 3562567..7f6d0c8 100644
--- a/polarization/tissueLayer.py
+++ b/polarization/tissueLayer.py
@@ -3,6 +3,7 @@
 from .pulse import Pulse
 from typing import List
 import numpy as np
+from collections import namedtuple
 
 __all__ = ['TissueLayer', 'SurfaceTissueLayer', 'RandomTissueLayer', 'EmptyTissueLayer']
 
@@ -25,7 +26,13 @@ def __init__(self, birefringence, opticAxis, scattDensity=0, thickness=200):
 
         self.cachedMatrices = {}
         self.apparentOpticAxis = None
-        self.scatterers = ScattererGroup(self.thickness, self.scattDensity)
+        self.scatterers = Scatterers(self.thickness, self.scattDensity)
+
+    def copy(self, thickness=None):
+        if thickness is None:
+            thickness = self.thickness
+        return TissueLayer(birefringence=self.birefringence, opticAxis=self.opticAxis,
+                           scattDensity=self.scattDensity, thickness=thickness)
 
     @property
     def opticAxis(self):
@@ -43,7 +50,7 @@ def opticAxis(self, vector):
     def orientation(self):
         """ Orientation in radians in X/Y coordinates from the Q/U plane of linear polarization. """
         if self.opticAxis[0] == 0:
-            return np.pi / 2 / 2
+            return np.pi / 4
         else:
             return np.arctan(self.opticAxis[1] / self.opticAxis[0]) / 2
 
@@ -90,10 +97,10 @@ def scatteringDeltaAt(self, K):
             dX, dY = np.zeros(len(K), dtype=complex), np.zeros(len(K), dtype=complex)
             K = np.asarray(K)
 
-        for scat in self.scatterers:
-            phi = 2 * scat.dz * K
-            dX += scat.strength * exp(1j * phi)
-            dY += scat.strength * exp(1j * phi * (1 + self.birefringence))
+        for dz, strength in zip(self.scatterers.dz, self.scatterers.strength):
+            phi = 2 * dz * K
+            dX += strength * exp(1j * phi)
+            dY += strength * exp(1j * phi * (1 + self.birefringence))
         return dX, dY
 
     def backscatteringMatrixAt(self, k):
@@ -109,39 +116,47 @@ def resetScatterers(self):
         return self.scatterers.reset()
 
 
-class Scatterer:
-    def __init__(self, max_dz):
-        self.dz = np.random.rand() * max_dz
-        self.strength = np.random.rand()
-        self.transferMatrix = None
+Scatterer = namedtuple('Scatterer', ['dz', 'strength'])
+
+
+class RandomScatterer:
+    def __new__(cls, max_dz):
+        return Scatterer(dz=np.random.rand() * max_dz, strength=np.random.rand())
 
 
-class ScattererGroup:
-    def __init__(self, length, density):
+class Scatterers:
+    def __init__(self, length, density=None, N=None):
         self.length = length
-        self.N = int(density * length)
+        if density:
+            self.N = int(density * length)
+        else:
+            self.N = N
 
-        self.scatterers = None
+        self.dz = []
+        self.strength = []
         self.reset()
 
     def reset(self):
-        self.scatterers = []
-        for _ in range(self.N):
-            self.scatterers.append(Scatterer(self.length))
-
-    @property
-    def dz(self) -> list:
-        return [scatterer.dz for scatterer in self.scatterers]
-
-    @property
-    def strength(self) -> list:
-        return [scatterer.strength for scatterer in self.scatterers]
+        self.dz = np.random.rand(self.N) * self.length
+        self.strength = np.random.rand(self.N)
 
     def __iter__(self):
-        return iter(self.scatterers)
+        self.n = 0
+        return self
+
+    def __next__(self) -> Scatterer:
+        if self.n < self.N:
+            scatterer = Scatterer(self.dz[self.n], self.strength[self.n])
+            self.n += 1
+            return scatterer
+        else:
+            raise StopIteration
 
     def __len__(self):
-        return len(self.scatterers)
+        return self.N
+
+    def __getitem__(self, item):
+        return Scatterer(self.dz[item], self.strength[item])
 
 
 class SurfaceTissueLayer(TissueLayer):
diff --git a/polarization/tissueStack.py b/polarization/tissueStack.py
index 6d80651..9b8f840 100644
--- a/polarization/tissueStack.py
+++ b/polarization/tissueStack.py
@@ -1,7 +1,7 @@
 from .tissueLayer import *
 from .jonesvector import JonesVector
 from .jonesmatrix import JonesMatrix, Vacuum
-from .pulse import Pulse
+from .pulse import *
 from typing import List
 import numpy as np
 
@@ -9,7 +9,7 @@
 
 
 class TissueStack:
-    def __init__(self, offset=0, layers=None):
+    def __init__(self, offset=0, layers=None, height=3000):
         """
         A stack of multiple tissue layers along depth axis.
 
@@ -18,11 +18,20 @@ def __init__(self, offset=0, layers=None):
         """
         self.layers: List[TissueLayer] = []
         self.offset = offset
+        self.height = height
+
+        self.forwardMatrices = None
+        self.backwardMatrices = None
 
         if layers is not None:
             for layer in layers:
                 self.append(layer)
 
+        self._scattDensity = None
+        self._opticAxis = None
+        self._apparentOpticAxis = None
+        self._birefringence = None
+
     def append(self, layer: TissueLayer):
         self.layers.append(layer)
 
@@ -32,19 +41,21 @@ def __len__(self):
     def __iter__(self):
         return iter(self.layers)
 
-    def transferMatrix(self, layerIndex=None, backward=False):
-        M = JonesMatrix(1, 0, 0, 1)
+    def initTransferMatrices(self):
+        self.forwardMatrices = [Vacuum(physicalLength=self.offset)]
+        self.backwardMatrices = [Vacuum(physicalLength=self.offset)]
+        for i, layer in enumerate(self.layers):
+            self.forwardMatrices.append(layer.transferMatrix() * self.forwardMatrices[i])
+            self.backwardMatrices.append(self.backwardMatrices[i] * layer.transferMatrix().backward())
+
+    def transferMatrix(self, layerIndex=-1, backward=False):
+        if self.forwardMatrices is None:
+            self.initTransferMatrices()
+
         if backward:
-            backwardLayers = self.layers[:layerIndex]
-            backwardLayers.reverse()
-            for layer in backwardLayers:
-                M *= layer.transferMatrix().backward()
-            M *= Vacuum(physicalLength=self.offset)
+            return self.backwardMatrices[layerIndex]
         else:
-            M *= Vacuum(physicalLength=self.offset)
-            for layer in self.layers[: layerIndex]:
-                M *= layer.transferMatrix()
-        return M
+            return self.forwardMatrices[layerIndex]
 
     def propagateThrough(self, vector: JonesVector) -> JonesVector:
         return self.transferMatrix() * vector
@@ -58,28 +69,94 @@ def propagateManyThrough(self, vectors: List[JonesVector]) -> List[JonesVector]:
     def backscatter(self, vector: JonesVector) -> JonesVector:
         signal = JonesVector(0, 0, k=vector.k)
         for i, layer in enumerate(self.layers):
-            signal += self.transferMatrix(i, backward=True) * (self.transferMatrix(i) * layer.backscatter(vector))
+            signal += self.transferMatrix(i, backward=True) * (layer.backscatter(self.transferMatrix(i) * vector))
         return signal
 
-    def backscatterMany(self, vectors: List[JonesVector]):
-        if type(vectors) is Pulse:
-            K = vectors.k
+    def backscatterMany(self, vectors):
+        if type(vectors) is PulseCollection:
+            return self._backscatterPulseCollection(vectors)
+        elif type(vectors) is Pulse:
+            return self._backscatterPulse(vectors)
         else:
-            K = [v.k for v in vectors]
+            return self._backscatterVectors(vectors)
+
+    def _backscatterVectors(self, vectors: List[JonesVector]):
+        K = [v.k for v in vectors]
         self.initBackscatteringAt(K)
 
         vectorsOut = []
         for v in vectors:
             vectorsOut.append(self.backscatter(v))
-
+        
         if type(vectors) is Pulse:
             return Pulse(vectors=vectorsOut)
         return vectorsOut
 
+    def _backscatterPulse(self, pulse: Pulse):
+        return self._backscatterPulseCollection(PulseCollection(pulses=[pulse]))[0]
+
+    def _backscatterPulseCollection(self, pulses: PulseCollection):
+        self.initBackscatteringAt(pulses.k)
+
+        self.initTransferMatrices()
+
+        vectorsOut = [[] for _ in pulses]
+        for i, k in enumerate(pulses.k):
+            vectorOut = [JonesVector(0, 0, k=k) for _ in pulses]
+            for j, layer in enumerate(self.layers):
+                M = self.transferMatrix(j, backward=True) * layer.backscatteringMatrixAt(k) * self.transferMatrix(j)
+                for p, pulse in enumerate(pulses):
+                    vectorOut[p] += M * pulse.vectors[i]
+            for p in range(len(pulses)):
+                vectorsOut[p].append(vectorOut[p])
+
+        pulsesOut = [Pulse(vectors) for vectors in vectorsOut]
+        return pulsesOut
+
     def initBackscatteringAt(self, K):
         for layer in self.layers:
             layer.initBackscatteringMatrixAt(K)
 
+    @property
+    def scattDensity(self):
+        if self._scattDensity is None:
+            self._scattDensity = np.zeros(int(self.height))
+            currentPosition = int(self.offset)
+            for layer in self.layers:
+                self._scattDensity[currentPosition: currentPosition + int(layer.thickness)] = layer.scattDensity
+                currentPosition += int(layer.thickness)
+        return self._scattDensity
+
+    @property
+    def opticAxis(self):
+        if self._opticAxis is None:
+            self._opticAxis = np.zeros((3, self.height))
+            currentPosition = int(self.offset)
+            for layer in self.layers:
+                self._opticAxis[:, currentPosition: currentPosition + int(layer.thickness)] = layer.opticAxis.reshape((3, 1))
+                currentPosition += int(layer.thickness)
+        return self._opticAxis
+
+    @property
+    def apparentOpticAxis(self):
+        if self._apparentOpticAxis is None:
+            self._apparentOpticAxis = np.zeros((3, self.height))
+            currentPosition = int(self.offset)
+            for layer in self.layers:
+                self._apparentOpticAxis[:, currentPosition: currentPosition + int(layer.thickness)] = layer.apparentOpticAxis.reshape((3, 1))
+                currentPosition += int(layer.thickness)
+        return self._apparentOpticAxis
+
+    @property
+    def birefringence(self):
+        if self._birefringence is None:
+            self._birefringence = np.zeros(int(self.height))
+            currentPosition = int(self.offset)
+            for layer in self.layers:
+                self._birefringence[currentPosition: currentPosition + int(layer.thickness)] = layer.birefringence
+                currentPosition += int(layer.thickness)
+        return self._birefringence
+
 
 class RandomTissueStack(TissueStack):
     def __init__(self, surface=True, maxBirefringence=0.0042, nLayers=None, offset=None, layerHeightRange=(60, 400)):
diff --git a/polarization/utils.py b/polarization/utils.py
index 1a9168b..4402502 100644
--- a/polarization/utils.py
+++ b/polarization/utils.py
@@ -113,3 +113,30 @@ def printModuleClasses(moduleName):
     for name, obj in inspect.getmembers(sys.modules[moduleName]):
         if inspect.isclass(obj) and obj.__module__.startswith(moduleName):
             print(obj)
+
+
+def createNoise(signal, SNR, mean=0):
+    """ White Gaussian Noise from a given SNR (not in dB).
+    If signal is complex, we double the SNR, create 2 independent normal distributions, and merge as complex.
+
+    Example:
+        >>> data = np.full(1000, 5 + 2j)
+        >>> data += createNoise(data, SNR=5)
+    """
+    assert SNR != 0
+
+    if np.iscomplexobj(signal):
+        SNR *= 2
+    power = np.abs(np.squeeze(signal)) ** 2
+    meanDecibel = 10 * np.log10(np.mean(power))
+    noiseMeanDecibel = meanDecibel - 10*np.log10(SNR)
+    noiseMeanPower = 10 ** (noiseMeanDecibel / 10)
+
+    if not np.iscomplexobj(signal):
+        noise = np.random.normal(mean, np.sqrt(noiseMeanPower), size=signal.shape)
+    else:
+        shape = list(signal.shape)
+        shape.append(2)
+        noise = np.random.normal(mean, np.sqrt(noiseMeanPower), size=tuple(shape)).view(np.complex128)
+        noise = np.reshape(noise, signal.shape)
+    return noise
diff --git a/polarization/vector.py b/polarization/vector.py
index eb7ba3d..d6f844a 100644
--- a/polarization/vector.py
+++ b/polarization/vector.py
@@ -3,11 +3,14 @@
 import sys
 
 class Vector:
-    def __init__(self, x: float = 0, y: float = 0, z: float = 0):
+    def __init__(self, x: float = 0, y: float = 0, z: float = 0, copy=None):
         """ Access properties with x,y and z.
         Internally for speed, read access done with self._x """
-
-        if isinstance(x, (int, float)):
+        if copy is not None:
+            self._x = copy.x
+            self._y = copy.y 
+            self._z = copy.z             
+        elif isinstance(x, (int, float)):
             self._x = x
             self._y = y 
             self._z = z
@@ -434,10 +437,10 @@ def normalizedDotProduct(self, vector):
 
 
 class ConstVector(Vector):
-    def __init__(self, x: float = 0, y: float = 0, z: float = 0):
+    def __init__(self, x: float = 0, y: float = 0, z: float = 0, copy=None):
         self._abs = 0
         self._norm = 0
-        Vector.__init__(self, x, y, z)
+        Vector.__init__(self, x, y, z, copy)
         self._abs = self.abs()
         self._norm = self.normSquared()
 


Test patch:
diff --git a/polarization/tests/comparisonTest.py b/polarization/tests/comparisonTest.py
new file mode 100644
index 0000000..82e2c13
--- /dev/null
+++ b/polarization/tests/comparisonTest.py
@@ -0,0 +1,294 @@
+import envtest
+from polarization.tissue import *
+from polarization.tissueStack import *
+from polarization.tissueLayer import *
+from polarization.pulse import *
+
+np.random.seed(521)
+
+
+def extractSV(S1, S2, S3):
+    SV1 = np.concatenate(
+        (S1[:, 0::2, np.newaxis], S2[:, 0::2, np.newaxis],
+         S3[:, 0::2, np.newaxis]), axis=2)
+    SV2 = np.concatenate(
+        (S1[:, 1::2, np.newaxis], S2[:, 1::2, np.newaxis],
+         S3[:, 1::2, np.newaxis]), axis=2)
+
+    return SV1, SV2
+
+
+def interweave(A, B, axis=0):
+    if axis != 0:
+        A = np.moveaxis(A, axis, 0)
+        B = np.moveaxis(B, axis, 0)
+
+    outShape = list(A.shape)
+    outShape[0] *= 2
+
+    C = np.empty(outShape, dtype=A.dtype)
+    C[0::2] = A
+    C[1::2] = B
+
+    if axis != 0:
+        C = np.moveaxis(C, 0, axis)
+    return C
+
+
+def euclideanNorm(arr, axis=0):
+    return np.expand_dims(np.linalg.norm(arr, axis=axis), axis=axis)
+
+
+def computeStokesFromFringes(fringes):
+    """ From fringes of shape (4, N*w, z) in k-space to Stokes vectors in z-space """
+    tomo = np.fft.fft(fringes, axis=-1)
+    return computeStokesFromTomo(tomo)
+
+
+def computeStokesFromTomo(tomo, weaveAxis=1):
+    """ From tomogram of shape (4, N*w, z) in z-space to Stokes vectors in z-space """
+    tomv1, tomh1, tomv2, tomh2 = tomo
+
+    tomoX = interweave(tomv1, tomv2, weaveAxis-1)  # (N*w*2, z)
+    tomoY = interweave(tomh1, tomh2, weaveAxis-1)
+
+    return computeStokesFromTomoXTomoY(tomoX, tomoY)
+
+
+def computeStokesFromTomoXTomoY(tomoX, tomoY):
+    S0 = np.absolute(tomoX) ** 2 + np.absolute(tomoY) ** 2
+    S1 = np.absolute(tomoX) ** 2 - np.absolute(tomoY) ** 2
+    S2 = 2 * np.real(tomoX * np.conj(tomoY))
+    S3 = -2 * np.imag(tomoX * np.conj(tomoY))
+
+    return np.asarray([S0, S1, S2, S3])
+
+
+class TestPSSignal(envtest.MyTestCase):
+    @staticmethod
+    def stokesFromFringes(fringes):
+        stokes = computeStokesFromFringes(fringes)
+        S0, S1, S2, S3 = np.moveaxis(stokes, 1, 2)
+
+        SV1, SV2 = extractSV(S1, S2, S3)
+        SV1 = SV1 / euclideanNorm(SV1, axis=2)
+        SV2 = SV2 / euclideanNorm(SV2, axis=2)
+
+        x = np.concatenate([SV1, SV2], axis=2)
+        return np.moveaxis(x, 0, 2)
+
+    @staticmethod
+    def createTestScan(layers):
+        resolution = 512
+        centerWavelength = 1.3
+        bandwidth = 0.13
+        Dk = 2 * np.pi / centerWavelength ** 2 * bandwidth
+        height = np.pi / Dk * resolution
+
+        tissue = TissueFromLayers(layers, height=int(height))
+        tissue.display()
+        # tissue.save("sim_eval_0_sample.npy", reSample=resolution)
+
+        pIn = PulseCollection.dualInputStates(centerWavelength, bandwidth, resolution=resolution)
+
+        return tissue.scan(pIn, verbose=True)
+
+    def runTestScan(self, layers):
+        pOut = self.createTestScan(layers)
+        # pOut.save("sim_eval_0_fringes.npy")
+
+        stokes = self.stokesFromFringes(pOut.fringes)
+
+        fig, axes = plt.subplots(1, 6)
+        for i, ax in enumerate(axes):
+            ax.imshow(stokes[:, i].T, aspect='auto', cmap='gray')
+        plt.show()
+
+    def testTissue(self):
+        """ Original Tissue configuration used to compare with old simulation.
+        Discontinuity in the interfaces around the 3rd layer (and following). """
+        layers = [TissueLayer(0.002, (0.1, 1, 0), 1000, 10), TissueLayer(0.0005, (1, -1, 0), 12, 600),
+                  TissueLayer(0.001, (0.2, 1, 0), 8, 500), TissueLayer(0.001, (-1, 0.3, 0), 3, 400),
+                  TissueLayer(0.0005, (0.8, -1, 0), 20, 500), TissueLayer(0.002, (1, 0.8, 0), 15, 300)]
+
+        # self.runTestScan(layers)
+
+    def testTissueUnit(self):
+        """ Smaller 4-layer tissue that highlights the same discontinuity. """
+        layers = [TissueLayer(0.0005, (1, -1, 0), 12, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.002, (-1, 0.3, 0), 3, 200),  # problematic layer
+                  TissueLayer(0.0005, (0.8, -1, 0), 20, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueSameLayers(self):
+        layers = [TissueLayer(0.001, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueScattering(self):
+        layers = [TissueLayer(0.001, (0.2, 1, 0), 12, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 3, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 20, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueBirefringence(self):
+        layers = [TissueLayer(0.0005, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.002, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.0005, (0.2, 1, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOA(self):
+        """ The discontinuity comes from OA changes. """
+        layers = [TissueLayer(0.001, (1, -1, 0), 8, 200),
+                  TissueLayer(0.001, (0.2, 1, 0), 8, 200),
+                  TissueLayer(0.001, (-1, 0.3, 0), 8, 200),
+                  TissueLayer(0.001, (0.8, -1, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOASimple(self):
+        layers = [TissueLayer(0.001, (1, -1, 0), 8, 200),
+                  TissueLayer(0.001, (0, 1, 0), 8, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (0, -1, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOAx(self):
+        layers = [TissueLayer(0.001, (1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOAy(self):
+        layers = [TissueLayer(0.001, (0, 1, 0), 8, 200),
+                  TissueLayer(0.001, (0, -1, 0), 8, 200),
+                  TissueLayer(0.001, (0, 1, 0), 8, 200),
+                  TissueLayer(0.001, (0, -1, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOASpin(self):
+        """ Rotate OA by 90 degrees each layer.
+        Classic discontinuity around 3rd layer. """
+        layers = [TissueLayer(0.001, (1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (0, 1, 0), 8, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (0, -1, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOASpin3Layers(self):
+        """ Discontinuity only before last layer. """
+        layers = [TissueLayer(0.001, (0, 1, 0), 8, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (0, -1, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOASpin3LayersThick(self):
+        """ Increasing thickness of 1st layer in a 3-layers configuration
+        still only creates a discontinuity before last layer. """
+        layers = [TissueLayer(0.001, (0, 1, 0), 1, 400),
+                  TissueLayer(0.001, (-1, 0, 0), 1, 200),
+                  TissueLayer(0.001, (0, -1, 0), 1, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOASpin2Layers(self):
+        """ No Discontinuity when only 2 layers are present. """
+        layers = [TissueLayer(0.001, (-1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (0, -1, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOAReverseSpin(self):
+        """ Reversed OA order doesn't change the problem. """
+        layers = [TissueLayer(0.001, (0, -1, 0), 8, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 200),
+                  TissueLayer(0.001, (0, 1, 0), 8, 200),
+                  TissueLayer(0.001, (1, 0, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testSmallTissueOASpin(self):
+        """ Discontinuity is visible at all interfaces after the 2nd layer. """
+        layers = [TissueLayer(0.001, (1, 0, 0), 8, 100),
+                  TissueLayer(0.001, (0, 1, 0), 8, 100),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 100),
+                  TissueLayer(0.001, (0, -1, 0), 8, 100),
+                  TissueLayer(0.001, (1, 0, 0), 8, 100),
+                  TissueLayer(0.001, (0, 1, 0), 8, 100),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 100),
+                  TissueLayer(0.001, (0, -1, 0), 8, 100)]
+
+        # self.runTestScan(layers)
+
+    def testSmallTissueOASlowSpin(self):
+        """ 45 degrees of rotation at each layer. """
+        layers = [TissueLayer(0.001, (1, 0, 0), 8, 100),
+                  TissueLayer(0.001, (1, 1, 0), 8, 100),
+                  TissueLayer(0.001, (0, 1, 0), 8, 100),
+                  TissueLayer(0.001, (-1, 1, 0), 8, 100),
+                  TissueLayer(0.001, (-1, 0, 0), 8, 100),
+                  TissueLayer(0.001, (-1, -1, 0), 8, 100),
+                  TissueLayer(0.001, (0, -1, 0), 8, 100)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOASlowerSpin(self):
+        """ 4 degrees or rotation at each layer. """
+        layers = [TissueLayer(0.001, (1, 0.0, 0), 8, 200),
+                  TissueLayer(0.001, (1, 0.1, 0), 8, 200),
+                  TissueLayer(0.001, (1, 0.2, 0), 8, 200),
+                  TissueLayer(0.001, (1, 0.3, 0), 8, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOASpinLowScat(self):
+        """ Discontinuity is of the same intensity (or same contrast at least) with few scatterers. """
+        layers = [TissueLayer(0.001, (1, 0, 0), 1, 200),
+                  TissueLayer(0.001, (0, 1, 0), 1, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 1, 200),
+                  TissueLayer(0.001, (0, -1, 0), 1, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueOASpinDecreaseFirstLayer(self):
+        """ The effect (around 3rd layer) is less visible when we decrease the size of the 1st layer. """
+        layers = [TissueLayer(0.001, (1, 0, 0), 1, 50),
+                  TissueLayer(0.001, (0, 1, 0), 1, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 1, 200),
+                  TissueLayer(0.001, (0, -1, 0), 1, 200)]
+
+        # self.runTestScan(layers)
+
+    def testTissueFastest(self):
+        """ Fastest configuration to see the effect. This is a copy of testTissueOASpinLowScat. """
+        layers = [TissueLayer(0.001, (1, 0, 0), 1, 200),
+                  TissueLayer(0.001, (0, 1, 0), 1, 200),
+                  TissueLayer(0.001, (-1, 0, 0), 1, 200),
+                  TissueLayer(0.001, (0, -1, 0), 1, 200)]
+
+        # self.runTestScan(layers)
+
+
+class TissueFromLayers(RandomTissue2D):
+    def __init__(self, layers, height):
+        testStack = TissueStack(offset=400, layers=layers)
+        super(TissueFromLayers, self).__init__(referenceStack=testStack, width=40, flat=True, height=height)
+
+
+if __name__ == '__main__':
+    envtest.main()
diff --git a/polarization/tests/testsJonesMatrix.py b/polarization/tests/testsJonesMatrix.py
index ce2aff6..d9bed24 100644
--- a/polarization/tests/testsJonesMatrix.py
+++ b/polarization/tests/testsJonesMatrix.py
@@ -5,25 +5,26 @@
 from numpy import exp, pi, angle, array, matmul, arctan2
 from numpy.linalg import eig, eigh
 
+
 class TestMatrices(envtest.MyTestCase):
     def testArrayRowCol(self):
-        m = array([[1,2],[3,4]])
+        m = array([[1, 2], [3, 4]])
         self.assertIsNotNone(m)
-        self.assertEqual(m[0,0], 1)
-        self.assertEqual(m[0,1], 2)
-        self.assertEqual(m[1,0], 3)
-        self.assertEqual(m[1,1], 4)
+        self.assertEqual(m[0, 0], 1)
+        self.assertEqual(m[0, 1], 2)
+        self.assertEqual(m[1, 0], 3)
+        self.assertEqual(m[1, 1], 4)
 
     def testMatrixProduct(self):
-        m = array([[1,2],[3,4]])
-        v = array([5,6])
-        r = matmul(m,v)
+        m = array([[1, 2], [3, 4]])
+        v = array([5, 6])
+        r = matmul(m, v)
         self.assertEqual(r[0], 17)
         self.assertEqual(r[1], 39)
 
     def testDefaultInitJonesMatrix(self):
         m = JonesMatrix()
-        
+
         self.assertIsNotNone(m)
         with self.assertRaises(ValueError) as context:
             m.m
@@ -34,129 +35,123 @@ def testDefaultInitJonesMatrix(self):
 
     def testDynamicPropertiesJonesMatrix(self):
         m = JonesMatrix()
-        
-        m.setValue('L', 2)        
-        self.assertEqual(m.L,2)
-        self.assertEqual(m.value('L'),2)
+
+        m.setValue('L', 2)
+        self.assertEqual(m.L, 2)
+        self.assertEqual(m.value('L'), 2)
 
     def testInitJonesMatrix(self):
-        m = JonesMatrix(1,0,0,1,physicalLength=1.0)
-        
+        m = JonesMatrix(1, 0, 0, 1, physicalLength=1.0)
+
         self.assertIsNotNone(m)
-        self.assertEqual(m.determinant,1)
+        self.assertEqual(m.determinant, 1)
         self.assertEqual(m.L, 1)
 
     def testInitJonesMatrixABCD(self):
-        m = JonesMatrix(1,2,3,4,physicalLength=1.0)
-        
+        m = JonesMatrix(1, 2, 3, 4, physicalLength=1.0)
+
         self.assertIsNotNone(m)
         self.assertEqual(m.L, 1)
+
         self.assertEqual(m.A, 1)
         self.assertEqual(m.B, 2)
         self.assertEqual(m.C, 3)
         self.assertEqual(m.D, 4)
 
     def testOrientedJonesMatrixABCD(self):
-        m = JonesMatrix(1,2,3,4,physicalLength=1.0)
-        determinant = det(m.m)
-        for theta in range(0,370,10):
-            m.orientation = theta*radPerDeg
-            self.assertAlmostEqual(det(m.m), determinant)
+        m = JonesMatrix(1, 2, 3, 4, physicalLength=1.0)
+        determinant = m.determinant
+        for theta in range(0, 370, 10):
+            m.orientation = theta * radPerDeg
+            self.assertAlmostEqual(m.determinant, determinant)
 
     def testFlippedJonesMatrixABCD(self):
-        m = JonesMatrix(1,2,3,4,physicalLength=1.0)
-        determinant = det(m.m)
-        m.orientation = 2*pi
-        self.assertAlmostEqual(det(m.m), determinant)
+        m = JonesMatrix(1, 2, 3, 4, physicalLength=1.0)
+        determinant = m.determinant
+        m.orientation = 2 * pi
+        self.assertAlmostEqual(m.determinant, determinant)
         self.assertAlmostEqual(m.A, 1)
         self.assertAlmostEqual(m.B, 2)
         self.assertAlmostEqual(m.C, 3)
         self.assertAlmostEqual(m.D, 4)
 
         m.orientation = pi
-        self.assertAlmostEqual(det(m.m), determinant)
+        self.assertAlmostEqual(m.determinant, determinant)
         self.assertAlmostEqual(m.A, 1)
         self.assertAlmostEqual(m.B, 2)
         self.assertAlmostEqual(m.C, 3)
         self.assertAlmostEqual(m.D, 4)
 
-    def testOrientedJonesMatrixABCD(self):
-        m = JonesMatrix(1,2,3,4,physicalLength=1.0)
-        determinant = det(m.m)
-        for theta in range(0,370,10):
-            m.orientation = theta*radPerDeg
-            self.assertAlmostEqual(det(m.m), determinant)
-
     def testPolarizerOrientation(self):
         h = HorizontalPolarizer()
         v = VerticalPolarizer()
         h.orientation = 0
-        v.orientation = -pi/2
+        v.orientation = -pi / 2
 
         p = LinearPolarizer(theta=0)
-        self.assertAlmostEqual(h.m[0,0], p.m[0,0])
-        self.assertAlmostEqual(h.m[0,1], p.m[0,1])
-        self.assertAlmostEqual(h.m[1,0], p.m[1,0])
-        self.assertAlmostEqual(h.m[1,1], p.m[1,1])
+        self.assertAlmostEqual(h.A, p.A)
+        self.assertAlmostEqual(h.B, p.B)
+        self.assertAlmostEqual(h.C, p.C)
+        self.assertAlmostEqual(h.D, p.D)
 
-        self.assertAlmostEqual(v.m[0,0], p.m[0,0])
-        self.assertAlmostEqual(v.m[0,1], p.m[0,1])
-        self.assertAlmostEqual(v.m[1,0], p.m[1,0])
-        self.assertAlmostEqual(v.m[1,1], p.m[1,1])
+        self.assertAlmostEqual(v.A, p.A)
+        self.assertAlmostEqual(v.B, p.B)
+        self.assertAlmostEqual(v.C, p.C)
+        self.assertAlmostEqual(v.D, p.D)
 
-        h.orientation = pi/2
+        h.orientation = pi / 2
         v.orientation = 0
 
-        p = LinearPolarizer(theta=pi/2)
-        self.assertAlmostEqual(h.m[0,0], p.m[0,0])
-        self.assertAlmostEqual(h.m[0,1], p.m[0,1])
-        self.assertAlmostEqual(h.m[1,0], p.m[1,0])
-        self.assertAlmostEqual(h.m[1,1], p.m[1,1])
+        p = LinearPolarizer(theta=pi / 2)
+        self.assertAlmostEqual(h.A, p.A)
+        self.assertAlmostEqual(h.B, p.B)
+        self.assertAlmostEqual(h.C, p.C)
+        self.assertAlmostEqual(h.D, p.D)
 
-        self.assertAlmostEqual(v.m[0,0], p.m[0,0])
-        self.assertAlmostEqual(v.m[0,1], p.m[0,1])
-        self.assertAlmostEqual(v.m[1,0], p.m[1,0])
-        self.assertAlmostEqual(v.m[1,1], p.m[1,1])
+        self.assertAlmostEqual(v.A, p.A)
+        self.assertAlmostEqual(v.B, p.B)
+        self.assertAlmostEqual(v.C, p.C)
+        self.assertAlmostEqual(v.D, p.D)
 
     def testMultiplyJonesMatrix(self):
-        m1 = JonesMatrix(1,0,0,1,physicalLength=1.0)
-        m2 = JonesMatrix(1,0,0,1,physicalLength=2.0)
-        
-        m = m1*m2
-        
+        m1 = JonesMatrix(1, 0, 0, 1, physicalLength=1.0)
+        m2 = JonesMatrix(1, 0, 0, 1, physicalLength=2.0)
+
+        m = m1 * m2
+
         self.assertIsNotNone(m)
-        self.assertEqual(m.determinant,1)
+        self.assertEqual(m.determinant, 1)
         self.assertEqual(m.L, 3)
 
     def testMultiplyJonesVectorForward(self):
-        m = JonesMatrix(1,2,3,4,physicalLength=1.0)
-        v = JonesVector(5,6)
+        m = JonesMatrix(1, 2, 3, 4, physicalLength=1.0)
+        v = JonesVector(5, 6)
         self.assertEqual(v.z, 0)
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertIsNotNone(vOut)
         self.assertEqual(vOut.z, 1.0)
 
     def testMultiplyJonesVectorBackward(self):
-        matrix = JonesMatrix(1,2,3,4,physicalLength=1.0)
-        v = JonesVector(5,6)
+        matrix = JonesMatrix(1, 2, 3, 4, physicalLength=1.0)
+        v = JonesVector(5, 6)
         self.assertEqual(v.z, 0)
 
         m = matrix.backward()
         v.reflect()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertIsNotNone(vOut)
         self.assertEqual(vOut.z, -1.0)
 
     def testTransformJonesVectorNoK(self):
-        m = JonesMatrix(1,2,3,4,physicalLength=1.0)
-        v = JonesVector(5,6)
+        m = JonesMatrix(1, 2, 3, 4, physicalLength=1.0)
+        v = JonesVector(5, 6)
         self.assertEqual(v.z, 0)
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertIsNotNone(vOut)
         self.assertEqual(vOut.Ex, 17)
@@ -165,11 +160,11 @@ def testTransformJonesVectorNoK(self):
         self.assertEqual(vOut.k, None)
 
     def testTransformJonesVectorWithK(self):
-        m = JonesMatrix(1,2,3,4,physicalLength=1.0)
-        v = JonesVector(5,6,k=6.28)
+        m = JonesMatrix(1, 2, 3, 4, physicalLength=1.0)
+        v = JonesVector(5, 6, k=6.28)
         self.assertEqual(v.z, 0)
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertIsNotNone(vOut)
         self.assertEqual(vOut.Ex, 17)
@@ -178,101 +173,101 @@ def testTransformJonesVectorWithK(self):
         self.assertEqual(vOut.k, 6.28)
 
     def testHorizontalPolarizer(self):
-        v = JonesVector(-1,1)
+        v = JonesVector(-1, 1)
         m = HorizontalPolarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertEqual(vOut.Ex, v.Ex)
         self.assertEqual(vOut.Ey, 0)
 
     def testHorizontalPolarizerBlock(self):
-        v = JonesVector(0,1)
+        v = JonesVector(0, 1)
         m = HorizontalPolarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertEqual(vOut.Ex, v.Ex)
         self.assertEqual(vOut.Ey, 0)
 
     def testVerticalPolarizer(self):
-        v = JonesVector(-1,1)
+        v = JonesVector(-1, 1)
         m = VerticalPolarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertEqual(vOut.Ex, 0)
         self.assertEqual(vOut.Ey, v.Ey)
 
     def testVerticalPolarizerBlock(self):
-        v = JonesVector(-1,0)
+        v = JonesVector(-1, 0)
         m = VerticalPolarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertEqual(vOut.Ex, 0)
         self.assertEqual(vOut.Ey, v.Ey)
 
     def testHorizontalPolarizerRotatedBy90GivesVertical(self):
-        h = HorizontalPolarizer().rotatedBy(pi/2)
+        h = HorizontalPolarizer().rotatedBy(pi / 2)
         v = VerticalPolarizer()
 
-        self.assertAlmostEqual(h.A,v.A)
-        self.assertAlmostEqual(h.B,v.B)
-        self.assertAlmostEqual(h.C,v.C)
-        self.assertAlmostEqual(h.D,v.D)
+        self.assertAlmostEqual(h.A, v.A)
+        self.assertAlmostEqual(h.B, v.B)
+        self.assertAlmostEqual(h.C, v.C)
+        self.assertAlmostEqual(h.D, v.D)
 
     def testPlus45Polarizer(self):
-        v = JonesVector(1,1)
+        v = JonesVector(1, 1)
         m = Plus45Polarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertEqual(vOut.Ex, v.Ex)
         self.assertEqual(vOut.Ey, v.Ey)
 
     def testPlus45PolarizerBlock(self):
-        v = JonesVector(-1,1)
+        v = JonesVector(-1, 1)
         m = Plus45Polarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertEqual(vOut.Ex, 0)
         self.assertEqual(vOut.Ey, 0)
 
     def testMinus45Polarizer(self):
-        v = JonesVector(-1,1)
+        v = JonesVector(-1, 1)
         m = Minus45Polarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertEqual(vOut.Ex, v.Ex)
         self.assertEqual(vOut.Ey, v.Ey)
 
     def testMinus45PolarizerBlock(self):
-        v = JonesVector(1,1)
+        v = JonesVector(1, 1)
         m = Minus45Polarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertEqual(vOut.Ex, 0)
         self.assertEqual(vOut.Ey, 0)
 
     def testRightCircularPolarizer(self):
-        v = JonesVector(1,1)
+        v = JonesVector(1, 1)
         m = RightCircularPolarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertTrue(vOut.isCircularlyPolarized)
         self.assertTrue(vOut.isRightCircularlyPolarized)
         self.assertFalse(vOut.isLeftCircularlyPolarized)
 
     def testLeftCircularPolarizer(self):
-        v = JonesVector(1,1)
+        v = JonesVector(1, 1)
         m = LeftCircularPolarizer()
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertTrue(vOut.isCircularlyPolarized)
         self.assertTrue(vOut.isLeftCircularlyPolarized)
@@ -286,14 +281,14 @@ def testRetarder(self):
         self.assertAlmostEqual(1, ret.D, 3)
 
     def testRetarderAsQWP(self):
-        ret = PhaseRetarder(delta=pi/2) # Ex is delayed by pi/2
-        self.assertAlmostEqual(exp(1j*pi/2), ret.A, 3)
+        ret = PhaseRetarder(delta=pi / 2)  # Ex is delayed by pi/2
+        self.assertAlmostEqual(exp(1j * pi / 2), ret.A, 3)
         self.assertAlmostEqual(0, ret.B, 3)
         self.assertAlmostEqual(0, ret.C, 3)
         self.assertAlmostEqual(1, ret.D, 3)
-        
-        v = JonesVector(1,1)
-        vOut = ret*v
+
+        v = JonesVector(1, 1)
+        vOut = ret * v
 
         self.assertTrue(vOut.isCircularlyPolarized)
         self.assertTrue(vOut.isRightCircularlyPolarized)
@@ -301,7 +296,7 @@ def testRetarderAsQWP(self):
 
     def testQuarterWaveplateInit(self):
         qwp = QWP(theta=0)
-        ret = PhaseRetarder(delta=-pi/2) # x is fast axis
+        ret = PhaseRetarder(delta=-pi / 2)  # x is fast axis
 
         self.assertAlmostEqual(qwp.A, ret.A, 3)
         self.assertAlmostEqual(qwp.B, ret.B, 3)
@@ -309,51 +304,51 @@ def testQuarterWaveplateInit(self):
         self.assertAlmostEqual(qwp.D, ret.D, 3)
 
     def testQuarterWaveplate(self):
-        v = JonesVector(1,1)
-        m = QWP(theta=0) # Ex is ahead = left circular
+        v = JonesVector(1, 1)
+        m = QWP(theta=0)  # Ex is ahead = left circular
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertTrue(vOut.isCircularlyPolarized)
         self.assertTrue(vOut.isLeftCircularlyPolarized)
         self.assertFalse(vOut.isRightCircularlyPolarized)
 
-        v = JonesVector(1,1)
-        m = QWP(theta=pi/2)
+        v = JonesVector(1, 1)
+        m = QWP(theta=pi / 2)
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertTrue(vOut.isCircularlyPolarized)
         self.assertTrue(vOut.isRightCircularlyPolarized)
         self.assertFalse(vOut.isLeftCircularlyPolarized)
 
     def testHalfWaveplate(self):
-        v = JonesVector(1,1)
-        m = HWP(theta=0) # Ex is ahead by pi
+        v = JonesVector(1, 1)
+        m = HWP(theta=0)  # Ex is ahead by pi
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertAlmostEqual(vOut.Ex.real, -1, 5)
         self.assertAlmostEqual(vOut.Ey.real, 1, 5)
 
-        v = JonesVector(1,1)
-        m = HWP(theta=pi/2) # Ey is ahead by pi
+        v = JonesVector(1, 1)
+        m = HWP(theta=pi / 2)  # Ey is ahead by pi
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertAlmostEqual(vOut.Ex.real, 1, 5)
         self.assertAlmostEqual(vOut.Ey.real, -1, 5)
 
-        v = JonesVector(1,0)
-        m = HWP(theta=pi/4) # +45 is ahead by pi/4
+        v = JonesVector(1, 0)
+        m = HWP(theta=pi / 4)  # +45 is ahead by pi/4
 
-        vOut = m*v
+        vOut = m * v
 
         self.assertAlmostEqual(vOut.Ex.real, 0, 5)
         self.assertAlmostEqual(vOut.Ey.real, -1, 5)
 
     def testRotation(self):
-        v = Rotation(theta=pi/2)*JonesVector(1,0)
+        v = Rotation(theta=pi / 2) * JonesVector(1, 0)
         self.assertAlmostEqual(abs(v.Ex), 0, 5)
         self.assertAlmostEqual(abs(v.Ey), 1, 5)
 
@@ -361,27 +356,27 @@ def testRotation(self):
     def testRetarders(self):
         phi, e1, e2 = PhaseRetarder(delta=0.4).birefringence
 
-        for theta in range(0,90,1):
-            M = PhaseRetarder(delta=0.4).rotatedBy(theta*radPerDeg)
+        for theta in range(0, 90, 1):
+            M = PhaseRetarder(delta=0.4).rotatedBy(theta * radPerDeg)
             delta, b1, b2 = M.birefringence
             if areAbsolutelyAlmostEqual(phi, delta):
                 # FIXME: There is a discontinuity at 45°
                 # It is not sufficient to flip b1 and b2
                 self.assertAlmostEqual(phi, delta)
-                self.assertAlmostEqual(b1.angleWith(e1, zHat), theta*radPerDeg)
+                self.assertAlmostEqual(b1.angleWith(e1, zHat), theta * radPerDeg)
             else:
                 self.assertTrue(False, "Discontinuity at: {0} {1} {2} {3}".format(theta, delta, b1, b2))
 
     def testBirefringenceInWaveplates(self):
         self.assertTrue(QWP(theta=0).isBirefringent)
-        self.assertTrue(QWP(theta=45*radPerDeg).isBirefringent)
-        self.assertTrue(QWP(theta=90*radPerDeg).isBirefringent)
-        self.assertTrue(QWP(theta=180*radPerDeg).isBirefringent)
+        self.assertTrue(QWP(theta=45 * radPerDeg).isBirefringent)
+        self.assertTrue(QWP(theta=90 * radPerDeg).isBirefringent)
+        self.assertTrue(QWP(theta=180 * radPerDeg).isBirefringent)
 
         self.assertTrue(HWP(theta=0).isBirefringent)
-        self.assertTrue(HWP(theta=45*radPerDeg).isBirefringent)
-        self.assertTrue(HWP(theta=90*radPerDeg).isBirefringent)
-        self.assertTrue(HWP(theta=180*radPerDeg).isBirefringent)
+        self.assertTrue(HWP(theta=45 * radPerDeg).isBirefringent)
+        self.assertTrue(HWP(theta=90 * radPerDeg).isBirefringent)
+        self.assertTrue(HWP(theta=180 * radPerDeg).isBirefringent)
 
     def testNoBirefringenceInOtherMatrices(self):
         self.assertFalse(HorizontalPolarizer().isBirefringent)
@@ -392,75 +387,76 @@ def testNoBirefringenceInOtherMatrices(self):
         self.assertFalse(LeftCircularPolarizer().isBirefringent)
 
     def testRotationNotBirefringent(self):
-        m = HorizontalPolarizer().rotatedBy(theta=pi/3)
+        m = HorizontalPolarizer().rotatedBy(theta=pi / 3)
         self.assertFalse(m.isBirefringent)
 
+    @unittest.skip("ProcketCell need to be checked")
     def testPockelsCell(self):
         c = PockelsCell(halfwaveVoltage=300, length=10)
         c.voltage = 150
-        self.assertAlmostEqual(c.retardance, pi/2)
+        self.assertAlmostEqual(c.retardance, pi / 2)
         c.voltage = -150
-        self.assertAlmostEqual(c.retardance, -pi/2)
+        self.assertAlmostEqual(c.retardance, -pi / 2)
 
     def testEigens(self):
-        for theta in [1,3,5,10,80,110,130,170,189]:
-            v1, v2, e1, e2 = Rotation(theta=theta*radPerDeg).eigens()
+        for theta in [1, 3, 5, 10, 80, 110, 130, 170, 189]:
+            v1, v2, e1, e2 = Rotation(theta=theta * radPerDeg).eigens()
             self.assertTrue(isNotZero(v1.imag), v1)
             self.assertTrue(isNotZero(v2.imag), v2)
-        for theta in [0,180]:
-            v1, v2, e1, e2 = Rotation(theta=theta*radPerDeg).eigens()
+        for theta in [0, 180]:
+            v1, v2, e1, e2 = Rotation(theta=theta * radPerDeg).eigens()
             self.assertTrue(isAlmostZero(v1.imag), v1)
             self.assertTrue(isAlmostZero(v2.imag), v2)
 
     def testArbitraryWavelengthDependentMatrix(self):
-        mat = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=0, physicalLength=1)
+        mat = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=0, physicalLength=1)
         with self.assertRaises(ValueError) as context:
             mat.m
 
         self.assertIsNotNone(mat.computeMatrix(k=2))
         anArray = mat.computeMatrix(k=2)
-        self.assertAlmostEqual( angle(anArray[1,1])-angle(anArray[0,0]), 0.2)
+        self.assertAlmostEqual(angle(anArray[1, 1]) - angle(anArray[0, 0]), 0.2)
 
     def testArbitraryWavelengthDependentMatrixAbrOrientation(self):
-        mat = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
+        mat = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
         with self.assertRaises(ValueError) as context:
             mat.m
-        self.assertEqual(mat.orientation, np.pi/2)
+        self.assertEqual(mat.orientation, np.pi / 2)
 
         self.assertIsNotNone(mat.computeMatrix(k=2))
         anArray = mat.computeMatrix(k=2)
-        self.assertAlmostEqual( angle(anArray[0,0])-angle(anArray[1,1]), 0.2)
+        self.assertAlmostEqual(angle(anArray[0, 0]) - angle(anArray[1, 1]), 0.2)
 
     def testBirefringentMaterialProductIsaMatrixProduct(self):
-        mat1 = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
-        mat2 = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
-        final = mat1*mat2
+        mat1 = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
+        mat2 = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
+        final = mat1 * mat2
         self.assertTrue(isinstance(final, MatrixProduct))
 
-        final2 = final*mat1
+        final2 = final * mat1
         self.assertTrue(isinstance(final2, MatrixProduct))
 
-        final3 = mat1*final
+        final3 = mat1 * final
         self.assertTrue(isinstance(final3, MatrixProduct))
 
     def testBirefringentMaterialProductOfMatrixProductIsAMatrixProduct(self):
-        mat1 = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
-        mat2 = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
-        final = mat1*mat2
-        final2 = final*mat1
+        mat1 = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
+        mat2 = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
+        final = mat1 * mat2
+        final2 = final * mat1
         self.assertTrue(isinstance(final2, MatrixProduct))
 
     def testMatrixProductProductWithAMatrixProductIsAMatrixProduct(self):
         mat1 = HorizontalPolarizer()
-        mat2 = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
-        mat3 = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
+        mat2 = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
+        mat3 = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
 
-        product12 = mat1*mat2
-        product23 = mat2*mat3
+        product12 = mat1 * mat2
+        product23 = mat2 * mat3
 
-        final = mat1*mat2*mat3
-        final2 = product12*mat3
-        final3 = mat1*product23
+        final = mat1 * mat2 * mat3
+        final2 = product12 * mat3
+        final3 = mat1 * product23
 
         self.assertEqual(final.matrices, final2.matrices)
         self.assertEqual(final.matrices, final3.matrices)
@@ -469,32 +465,32 @@ def testMatrixProductProductWithAMatrixProductIsAMatrixProduct(self):
         self.assertEqual(final.matrices[2], mat1)
 
     def testBirefringentMaterialProductVector(self):
-        mat = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
-        v = JonesVector(1,1,k=6.28)
-        vOut = mat*v
+        mat = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
+        v = JonesVector(1, 1, k=6.28)
+        vOut = mat * v
         self.assertIsNotNone(vOut)
         self.assertEqual(v.k, vOut.k)
 
     def testBirefringentMaterialMatrixProductVector(self):
         mat1 = HorizontalPolarizer()
-        mat2 = BirefringentMaterial(deltaIndex = 0.1, fastAxisOrientation=np.pi/2, physicalLength=1)
-        mat3 = BirefringentMaterial(deltaIndex = 0.2, fastAxisOrientation=np.pi/2, physicalLength=1)
-        v = JonesVector(1,0,k=2)
-        vOut = mat3*v
-        self.assertTrue(angle(vOut.Ex), 2*0.2)
+        mat2 = BirefringentMaterial(deltaIndex=0.1, fastAxisOrientation=np.pi / 2, physicalLength=1)
+        mat3 = BirefringentMaterial(deltaIndex=0.2, fastAxisOrientation=np.pi / 2, physicalLength=1)
+        v = JonesVector(1, 0, k=2)
+        vOut = mat3 * v
+        self.assertTrue(angle(vOut.Ex), 2 * 0.2)
         self.assertEqual(v.k, vOut.k)
-        vOut = mat2*mat3*v
-        self.assertTrue(angle(vOut.Ex), 2*0.3)
+        vOut = mat2 * mat3 * v
+        self.assertTrue(angle(vOut.Ex), 2 * 0.3)
         self.assertEqual(v.k, vOut.k)
-        product = mat2*mat3
-        vOut = product*v
-        self.assertTrue(angle(vOut.Ex), 2*0.3)
+        product = mat2 * mat3
+        vOut = product * v
+        self.assertTrue(angle(vOut.Ex), 2 * 0.3)
         self.assertEqual(v.k, vOut.k)
 
-
     def testBackwardMatrix(self):
         mat1 = HorizontalPolarizer()
         self.assertIsNotNone(mat1.backward())
 
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/polarization/tests/testsJonesVector.py b/polarization/tests/testsJonesVector.py
index 023903d..6f019e5 100644
--- a/polarization/tests/testsJonesVector.py
+++ b/polarization/tests/testsJonesVector.py
@@ -6,95 +6,96 @@
 import matplotlib
 import matplotlib.pyplot as plt
 
+
 class TestVector(envtest.MyTestCase):
     def testInitJones(self):
         v = JonesVector()
         self.assertIsNotNone(v)
 
-        v = JonesVector(1,0)
+        v = JonesVector(1, 0)
         self.assertIsNotNone(v)
 
-        v = JonesVector(complex(1,0),0)
+        v = JonesVector(complex(1, 0), 0)
         self.assertIsNotNone(v)
 
-        v = JonesVector(0, complex(1,0))
+        v = JonesVector(0, complex(1, 0))
         self.assertIsNotNone(v)
 
-        v = JonesVector(complex(1,0), complex(1,0))
+        v = JonesVector(complex(1, 0), complex(1, 0))
         self.assertIsNotNone(v)
-        self.assertEqual(v.z,0)
+        self.assertEqual(v.z, 0)
 
     def testComplexFunctions(self):
         n = complex(exp(-1j))
         self.assertEqual(angle(n), -1)
         n = complex(exp(10j))
-        self.assertAlmostEqual(angle(n), 10 - pi*4,3)
+        self.assertAlmostEqual(angle(n), 10 - pi * 4, 3)
 
     def testComponentsAreComplex(self):
-        v = JonesVector(complex(1,0), complex(1,0))
+        v = JonesVector(complex(1, 0), complex(1, 0))
         self.assertTrue(isinstance(v.Ex, complex))
         self.assertTrue(isinstance(v.Ey, complex))
 
-        v = JonesVector(1,1)
+        v = JonesVector(1, 1)
         self.assertTrue(isinstance(v.Ex, complex))
         self.assertTrue(isinstance(v.Ey, complex))
 
     def testStokesComponentsIntensity(self):
         v = JonesVector(1, 0)
-        self.assertEqual(v.S0,1)
+        self.assertEqual(v.S0, 1)
 
-        v = JonesVector(0,1)
-        self.assertEqual(v.S0,1)
+        v = JonesVector(0, 1)
+        self.assertEqual(v.S0, 1)
 
-        v = JonesVector(1,1)
-        self.assertEqual(v.S0,2)
+        v = JonesVector(1, 1)
+        self.assertEqual(v.S0, 2)
 
-        v = JonesVector(exp(1j),exp(2j))
-        self.assertEqual(v.S0,2)
+        v = JonesVector(exp(1j), exp(2j))
+        self.assertEqual(v.S0, 2)
 
-        v = JonesVector(exp(1j),exp(-2j))
-        self.assertEqual(v.S0,2)
+        v = JonesVector(exp(1j), exp(-2j))
+        self.assertEqual(v.S0, 2)
 
-        v = JonesVector(exp(-1j),exp(2j))
-        self.assertEqual(v.S0,2)
+        v = JonesVector(exp(-1j), exp(2j))
+        self.assertEqual(v.S0, 2)
 
-        v = JonesVector(-exp(-1j),-exp(2j))
-        self.assertEqual(v.S0,2)
+        v = JonesVector(-exp(-1j), -exp(2j))
+        self.assertEqual(v.S0, 2)
 
     def testIntensity(self):
         v = JonesVector(1, 1)
-        self.assertEqual(v.intensity,2)
+        self.assertEqual(v.intensity, 2)
 
     def testStokesComponentS0S1S2(self):
         v = JonesVector(1, 0)
-        self.assertEqual(v.S0,1)
-        self.assertEqual(v.S1,1)
-        self.assertEqual(v.S2,0)
+        self.assertEqual(v.S0, 1)
+        self.assertEqual(v.S1, 1)
+        self.assertEqual(v.S2, 0)
 
         v = JonesVector(0, 1)
-        self.assertEqual(v.S0,1)
-        self.assertEqual(v.S1,-1)
-        self.assertEqual(v.S2,0)
+        self.assertEqual(v.S0, 1)
+        self.assertEqual(v.S1, -1)
+        self.assertEqual(v.S2, 0)
 
         v = JonesVector(1, 1)
-        self.assertEqual(v.S0,2)
-        self.assertEqual(v.S1,0)
-        self.assertEqual(v.S2,2)
+        self.assertEqual(v.S0, 2)
+        self.assertEqual(v.S1, 0)
+        self.assertEqual(v.S2, 2)
 
         v = JonesVector(-1, -1)
-        self.assertEqual(v.S0,2)
-        self.assertEqual(v.S1,0)
-        self.assertEqual(v.S2,2)
+        self.assertEqual(v.S0, 2)
+        self.assertEqual(v.S1, 0)
+        self.assertEqual(v.S2, 2)
 
         v = JonesVector(1, -1)
-        self.assertEqual(v.S0,2)
-        self.assertEqual(v.S1,0)
-        self.assertEqual(v.S2,-2)
+        self.assertEqual(v.S0, 2)
+        self.assertEqual(v.S1, 0)
+        self.assertEqual(v.S2, -2)
 
         v = JonesVector(-1, 1)
-        self.assertEqual(v.S0,2)
-        self.assertEqual(v.S1,0)
-        self.assertEqual(v.S2,-2)
+        self.assertEqual(v.S0, 2)
+        self.assertEqual(v.S1, 0)
+        self.assertEqual(v.S2, -2)
 
     def testNormalization(self):
         v = JonesVector(1, 1)
@@ -116,7 +117,6 @@ def testStokesComponentS0S1S2AreReal(self):
         with self.assertRaises(NotImplementedError):
             v.StokesVector
 
-
     def testLinearPolarization(self):
         v = JonesVector(1, 0)
         self.assertTrue(v.isLinearlyPolarized)
@@ -160,36 +160,36 @@ def testCircularPolarization(self):
         v = JonesVector(1, 0)
         self.assertFalse(v.isCircularlyPolarized)
 
-        v = JonesVector(1, exp(-1j*pi/2))
+        v = JonesVector(1, exp(-1j * pi / 2))
         self.assertTrue(v.isRightCircularlyPolarized)
 
-        v = JonesVector(1, exp(1j*pi/2))
+        v = JonesVector(1, exp(1j * pi / 2))
         self.assertTrue(v.isLeftCircularlyPolarized)
 
-        v = JonesVector(exp(1j*pi/2), 1)
+        v = JonesVector(exp(1j * pi / 2), 1)
         self.assertTrue(v.isRightCircularlyPolarized)
 
-        v = JonesVector(exp(-1j*pi/2), 1)
+        v = JonesVector(exp(-1j * pi / 2), 1)
         self.assertTrue(v.isLeftCircularlyPolarized)
 
-        v = JonesVector(exp(-1j*pi/2), 0.9)
+        v = JonesVector(exp(-1j * pi / 2), 0.9)
         self.assertFalse(v.isCircularlyPolarized)
         self.assertTrue(v.isEllipticallyPolarized)
 
     def testCircularPolarizationLargePhase(self):
-        v = JonesVector(1, exp(-1j*5*pi/2))
+        v = JonesVector(1, exp(-1j * 5 * pi / 2))
         self.assertTrue(v.isRightCircularlyPolarized)
 
-        v = JonesVector(1, exp(1j*5*pi/2))
+        v = JonesVector(1, exp(1j * 5 * pi / 2))
         self.assertTrue(v.isLeftCircularlyPolarized)
 
-        v = JonesVector(exp(1j*5*pi/2), 1)
+        v = JonesVector(exp(1j * 5 * pi / 2), 1)
         self.assertTrue(v.isRightCircularlyPolarized)
 
-        v = JonesVector(exp(-1j*5*pi/2), 1)
+        v = JonesVector(exp(-1j * 5 * pi / 2), 1)
         self.assertTrue(v.isLeftCircularlyPolarized)
 
-        v = JonesVector(exp(-1j*5*pi/2), 0.9)
+        v = JonesVector(exp(-1j * 5 * pi / 2), 0.9)
         self.assertFalse(v.isCircularlyPolarized)
         self.assertTrue(v.isEllipticallyPolarized)
 
@@ -227,37 +227,37 @@ def testPredefined(self):
         self.assertTrue(v.isLinearlyPolarized)
 
     def testOrientation(self):
-        self.assertAlmostEqual(JonesVector.vertical().orientation*degPerRad, 90)
-        self.assertAlmostEqual(JonesVector.horizontal().orientation*degPerRad, 0)
-        self.assertAlmostEqual(JonesVector.plus45().orientation*degPerRad, 45)
-        self.assertAlmostEqual(JonesVector.minus45().orientation*degPerRad, -45)
-        self.assertAlmostEqual(JonesVector.rightCircular().orientation*degPerRad, 45)
-        self.assertAlmostEqual(JonesVector.leftCircular().orientation*degPerRad, 45)
+        self.assertAlmostEqual(JonesVector.vertical().orientation * degPerRad, 90)
+        self.assertAlmostEqual(JonesVector.horizontal().orientation * degPerRad, 0)
+        self.assertAlmostEqual(JonesVector.plus45().orientation * degPerRad, 45)
+        self.assertAlmostEqual(JonesVector.minus45().orientation * degPerRad, -45)
+        self.assertAlmostEqual(JonesVector.rightCircular().orientation * degPerRad, 45)
+        self.assertAlmostEqual(JonesVector.leftCircular().orientation * degPerRad, 45)
 
     def testGetSetValue(self):
         v = JonesVector(1, 1)
-        self.assertEqual(v.value('intensity'),2)
-        v.setValue('E1',0)
-        self.assertEqual(v.value('intensity'),1)
+        self.assertEqual(v.value('intensity'), 2)
+        v.setValue('E1', 0)
+        self.assertEqual(v.value('intensity'), 1)
 
     def testPhysicalField(self):
         v = JonesVector(1, 0)
         self.assertAlmostEqual(v.physicalField(0)[0], 1)
         self.assertAlmostEqual(v.physicalField(0)[1], 0)
-        self.assertAlmostEqual(v.physicalField(pi/2)[0], 0)
-        self.assertAlmostEqual(v.physicalField(pi/2)[1], 0)
+        self.assertAlmostEqual(v.physicalField(pi / 2)[0], 0)
+        self.assertAlmostEqual(v.physicalField(pi / 2)[1], 0)
 
         v = JonesVector.rightCircular()
-        self.assertAlmostEqual(v.physicalField(0)[0], 1/sqrt(2))
+        self.assertAlmostEqual(v.physicalField(0)[0], 1 / sqrt(2))
         self.assertAlmostEqual(v.physicalField(0)[1], 0)
-        self.assertAlmostEqual(v.physicalField(pi/2)[0], 0)
-        self.assertAlmostEqual(v.physicalField(pi/2)[1], 1/sqrt(2))
+        self.assertAlmostEqual(v.physicalField(pi / 2)[0], 0)
+        self.assertAlmostEqual(v.physicalField(pi / 2)[1], 1 / sqrt(2))
 
         v = JonesVector.leftCircular()
-        self.assertAlmostEqual(v.physicalField(0)[0], 1/sqrt(2))
+        self.assertAlmostEqual(v.physicalField(0)[0], 1 / sqrt(2))
         self.assertAlmostEqual(v.physicalField(0)[1], 0)
-        self.assertAlmostEqual(v.physicalField(pi/2)[0], 0)
-        self.assertAlmostEqual(v.physicalField(pi/2)[1], -1/sqrt(2))
+        self.assertAlmostEqual(v.physicalField(pi / 2)[0], 0)
+        self.assertAlmostEqual(v.physicalField(pi / 2)[1], -1 / sqrt(2))
 
     def testPhysicalRealField(self):
         v = JonesVector(1, 0)
@@ -269,11 +269,11 @@ def testFormat(self):
         self.assertEqual("Ex = 1.00, Ey = 0.00", "{0}".format(v))
         v = JonesVector(0, 1)
         self.assertEqual("Ex = 0.00, Ey = 1.00", "{0}".format(v))
-        v = JonesVector(exp(1j*pi/2), 1)
+        v = JonesVector(exp(1j * pi / 2), 1)
         self.assertEqual("Ex = exp(π/2j), Ey = 1.00", "{0}".format(v))
-        v = JonesVector(1, exp(-1j*pi/2))
+        v = JonesVector(1, exp(-1j * pi / 2))
         self.assertEqual("Ex = 1.00, Ey = exp(-π/2j)", "{0}".format(v))
-        v = JonesVector(0.5*exp(-1j*pi/2), 0.1*exp(1j*pi/4))
+        v = JonesVector(0.5 * exp(-1j * pi / 2), 0.1 * exp(1j * pi / 4))
         self.assertEqual("Ex = 0.50 ⨉ exp(-π/2j), Ey = 0.10 ⨉ exp(π/4j)", "{0}".format(v))
 
     def testReflection(self):
@@ -284,5 +284,6 @@ def testReflection(self):
         self.assertTrue(v.isLeftCircularlyPolarized)
         self.assertEqual(v.b1.cross(v.b2), v.b3)
 
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/polarization/tests/testsTissue.py b/polarization/tests/testsTissue.py
index d1daa77..2e9f95d 100644
--- a/polarization/tests/testsTissue.py
+++ b/polarization/tests/testsTissue.py
@@ -12,8 +12,8 @@ def setUp(self) -> None:
         pass
 
     def testRandomTissue2DMap(self):
-        tissue = RandomTissue2D(nLayers=6, width=5)
-        self.assertEqual(tissue._layerSizeMap.shape, (7, 5))
+        tissue = RandomTissue2D(nLayers=6)
+        self.assertEqual(tissue._layerSizeMap.shape, (7, 200))
 
     def testPSOCTFringes(self):
         resolution = 5
@@ -44,7 +44,8 @@ def testPSOCT(self):
 
 class TissueTestUnit(RandomTissue2D):
     def __init__(self):
-        layers = [TissueLayer(0.004, (0, 1, 0), 20, 300), TissueLayer(0.004, (1, 0, 0), 1, 200)]
+        layers = [TissueLayer(0.004, (0, 1, 0), 20, 300), TissueLayer(0.004, (1, 0, 0), 1, 200),
+                  TissueLayer(0.001, (0, 1, 0), 20, 300), TissueLayer(0.001, (1, 0, 0), 1, 200)]
         testStack = TissueStack(offset=300, layers=layers)
         super(TissueTestUnit, self).__init__(referenceStack=testStack, width=4, flat=True)
 
diff --git a/polarization/tests/testsTissueLayer.py b/polarization/tests/testsTissueLayer.py
index 79ae2bc..0c68f76 100644
--- a/polarization/tests/testsTissueLayer.py
+++ b/polarization/tests/testsTissueLayer.py
@@ -1,6 +1,6 @@
 import envtest
 from polarization import *
-from polarization.tissueLayer import Scatterer, ScattererGroup
+from polarization.tissueLayer import Scatterers
 
 np.random.seed(521)
 
@@ -30,8 +30,8 @@ def testSymetricTransferMatrix(self):
 
         MRef = self.layerRef.transferMatrix(k=self.k)
 
-        for a, b in zip([J.A, J.B, J.C, J.D], np.nditer(MRef)):
-            self.assertAlmostEqual(a, b, 10)
+        for a, b in zip(J.m, np.nditer(MRef)):
+            self.assertAlmostEqual(a, b)
 
     def testPropagationTransferMatrixNotEqual(self):
         """ Not equal because we are not defining our retarder matrix as symmetric. """
@@ -40,7 +40,7 @@ def testPropagationTransferMatrixNotEqual(self):
         MRef = self.layerRef.transferMatrix(k=self.k)
 
         for a, b in zip(np.nditer(M), np.nditer(MRef)):
-            self.assertNotAlmostEqual(a, b, 10)
+            self.assertNotAlmostEqual(a, b)
 
     def testPropagateWithPhaseSymetricMatrix(self):
         delta = self.k * self.birefringence * self.thickness
@@ -51,9 +51,12 @@ def testPropagateWithPhaseSymetricMatrix(self):
 
         MRef = self.layerRef.transferMatrix(k=self.k)
         pOutRef = np.reshape(np.einsum('ij, j', MRef, np.asarray([self.pIn.Ex, self.pIn.Ey])), (2,))
+        pOutRef = JonesVector(pOutRef[0], pOutRef[1])
+
+        self.assertAlmostEqual(pOut.Ex, pOutRef.Ex)
+        self.assertAlmostEqual(pOut.Ey, pOutRef.Ey)
 
-        self.assertAlmostEqual(pOut.Ex, pOutRef[0])
-        self.assertAlmostEqual(pOut.Ey, pOutRef[1])
+        self.assertAlmostEqual(pOut.orientation, pOutRef.orientation)
 
     def testPropagate(self):
         """ When using our not phase-symmetric matrix, we expect same output orientation, but different phase."""
@@ -63,9 +66,6 @@ def testPropagate(self):
         pOutRef = np.reshape(np.einsum('ij, j', MRef, np.asarray([self.pIn.Ex, self.pIn.Ey])), (2,))
         pOutRef = JonesVector(pOutRef[0], pOutRef[1])
 
-        self.assertNotAlmostEqual(pOut.Ex, pOutRef.Ex)
-        self.assertNotAlmostEqual(pOut.Ey, pOutRef.Ey)
-
         self.assertAlmostEqual(pOut.orientation, pOutRef.orientation)
 
     def testPropagateMany(self):
@@ -83,7 +83,7 @@ def testPropagateBackward(self):
         pOutRef = np.reshape(np.einsum('ij, j', MRef, np.asarray([self.pIn.Ex, self.pIn.Ey])), (2,))
         pOutRef = JonesVector(pOutRef[0], pOutRef[1])
 
-        self.assertAlmostEqual(pOut.orientation, pOutRef.orientation, places=10)
+        self.assertAlmostEqual(pOut.orientation, pOutRef.orientation)
 
     def testScatterers(self):
         scatStrength = self.layer.scatterers.strength
@@ -97,7 +97,7 @@ def testScatterers(self):
         self.assertTrue(np.max(scatDz) <= self.thickness)
 
     def testSingleBackscatter(self):
-        self.layer.scatterers = [Scatterer(self.thickness)]
+        self.layer.scatterers = Scatterers(self.thickness, N=1)
 
         pOut = self.layer.backscatter(self.pIn)
 
@@ -113,7 +113,7 @@ def testBackscatterSum(self):
         """ Summing the backscattered signal of 2 scatterers. While each scatSignal and reference signals have
         the same orientations, their sum doesn't because the phases are different. """
         # a phase-symmetric matrix gives the same results. Test with backward() instead of 2*dz when possible
-        self.layer.scatterers = [Scatterer(self.thickness), Scatterer(self.thickness)]
+        self.layer.scatterers = Scatterers(self.thickness, N=2)
 
         # TissueLayer
         pOut = self.layer.backscatter(self.pIn)
diff --git a/polarization/tests/testsTissueStack.py b/polarization/tests/testsTissueStack.py
index 72d2e57..2455191 100644
--- a/polarization/tests/testsTissueStack.py
+++ b/polarization/tests/testsTissueStack.py
@@ -1,56 +1,264 @@
 import envtest
 from polarization import *
+from polarization.tissueLayer import Scatterers
 
 np.random.seed(521)
 
 
-class TestTissueStack(envtest.MyTestCase):
-    def setUp(self) -> None:
-        self.stack = TissueStackTestUnit()
+class TestTissuePropagation(envtest.MyTestCase):
+    def __init__(self, *args, **kwargs):
+        self.k = None
+        self.pIn = None
+        super(TestTissuePropagation, self).__init__(*args, **kwargs)
 
+    def setUp(self):
         self.k = 2 * np.pi / 1.3
         self.pIn = JonesVector.horizontal()
         self.pIn.k = self.k
 
-    def testPropagate(self):
-        """ When using our not phase-symmetric matrix, we expect same output orientation, but different phase."""
-        pOut = self.stack.propagateThrough(self.pIn)
+    def testPropagateThroughNoBirefringence(self):
+        """ Should equal propagation in Vacuum. """
+        stack = TissueStackNoBirefringence()
+        stackLength = stack.offset + sum([layer.thickness for layer in stack.layers])
 
-        print(pOut.orientation)
+        pOut = stack.propagateThrough(self.pIn)
 
-        # no ref to compare with
+        self.assertAlmostEqual(pOut.Ex, exp(1j * self.pIn.k * stackLength))
+        self.assertAlmostEqual(pOut.Ey, 0)
 
-    def testPropagateMany(self):
+    def testPropagateThroughOrientedTissueHorizontal(self):
+        """ With all layer optic axes oriented with the beam (no Q/U components), there should be no retarding effect
+        other than Vacuum propagation. """
+        stack = TissueStackOrientedHorizontally()
+        stackLength = stack.offset + sum([layer.thickness for layer in stack.layers])
+
+        pOut = stack.propagateThrough(self.pIn)
+
+        self.assertAlmostEqual(pOut.Ex, exp(1j * self.pIn.k * stackLength))
+        self.assertAlmostEqual(pOut.Ey, 0)
+
+    def testPropagateThroughOrientedTissue45(self):
+        """ With all layer optic axes oriented with the beam (no Q/U components), there should be no retarding effect
+        other than Vacuum propagation. """
+        self.pIn = JonesVector.plus45()
+        self.pIn.k = self.k
+        stack = TissueStackOriented45Degrees()
+        stackLength = stack.offset + sum([layer.thickness for layer in stack.layers])
+
+        pOut = stack.propagateThrough(self.pIn)
+
+        self.assertAlmostEqual(pOut.Ex, sin(np.pi / 4) * exp(1j * self.pIn.k * stackLength))
+        self.assertAlmostEqual(pOut.Ey, sin(np.pi / 4) * exp(1j * self.pIn.k * stackLength))
+
+    def testPropagateThroughPerpendicularLayer(self):
+        """ With a beam polarization perpendicular to the one and only tissue optic axis,
+        the retarding effect is simply ikL(1+dn). """
+        self.pIn = JonesVector.vertical()
+        self.pIn.k = self.k
+        stack = TissueStackSingleHorizontal()
+        stackLength = stack.offset + sum([layer.thickness for layer in stack.layers])
+        layerLength = stack.layers[0].thickness
+        layerBirefringence = stack.layers[0].birefringence
+
+        pOut = stack.propagateThrough(self.pIn)
+
+        self.assertAlmostEqual(pOut.Ex, 0)
+        self.assertAlmostEqual(pOut.Ey, exp(1j * self.pIn.k * (stackLength + layerLength * layerBirefringence)))
+
+    def testPropagateThroughSingleLayer(self):
+        """ With a single arbitrary birefringent layer, the test becomes more involved.
+        For a single rotation matrix on vacuum propagation, we have [[A, B], [C, D]] * [e^ikl, 0]^T
+        """
+        stack = TissueStackSingle()
+        layerLength = stack.layers[0].thickness
+        layerBirefringence = stack.layers[0].birefringence
+
+        jPhi = 1j * layerLength * self.pIn.k
+        A = exp(jPhi)
+        D = exp(jPhi * (1 + layerBirefringence))
+        orientation = stack.layers[0].orientation
+        Ap = cos(orientation) ** 2 * A + D * sin(orientation) ** 2
+        Bp = (A - D) * sin(orientation) * cos(orientation)
+        Dp = cos(orientation) ** 2 * D + A * sin(orientation) ** 2
+
+        pOut = stack.propagateThrough(self.pIn)
+
+        self.assertAlmostEqual(pOut.Ex, Ap * exp(1j * self.pIn.k * stack.offset))
+        self.assertAlmostEqual(pOut.Ey, Bp * exp(1j * self.pIn.k * stack.offset))
+
+    def testPropagateThrough(self):
+        """ With multiple arbitrary birefringent layers, we have to use the results from last layer
+        and since the input at secondary layers is not horizontal, the math is a little more complex.
+        """
+        # todo, maybe with numpy.matmul ...
+        # fisrt layer : [[A, B], [C, D]] * [e^ikl, 0]^T = e^ikL * [A , B]^T
+        # second layer: [[E, F], [G, H]] * e^ikL * [A, B]^T = ...
+        pass
+
+    def testPropagateManyThrough(self):
+        """ Qualitative test for pulse propagation (multiple vectors). """
+        # todo: quantitative test
+
+        stack = TissueStackUnit()
         res = 5
         pIn = Pulse.horizontal(centerWavelength=1.3, wavelengthBandwidth=0.13, resolution=res)
-        pOut = self.stack.propagateManyThrough(pIn)
+        pOut = stack.propagateManyThrough(pIn)
 
         self.assertTrue(len(pOut) == res)
-        self.assertTrue(pOut[0].orientation != pOut[res//2].orientation)
+        self.assertTrue(pOut[0].orientation != pOut[res // 2].orientation)
 
-    def testBackscatter(self):
-        pOut = self.stack.backscatter(self.pIn)
 
-        print(pOut.orientation)
+class TestTissueBackscattering(envtest.MyTestCase):
+    def __init__(self, *args, **kwargs):
+        self.k = None
+        self.pIn = None
+        super(TestTissueBackscattering, self).__init__(*args, **kwargs)
+
+    def setUp(self):
+        self.k = 2 * np.pi / 1.3
+        self.pIn = JonesVector.horizontal()
+        self.pIn.k = self.k
+
+    def testBackscatterSingleScattererSingleLayerNoBirefringence(self):
+        stack = TissueStackSingleScattererSingleLayerNoBirefringence()
+
+        scatPosition = stack.offset + stack.layers[0].scatterers.dz[0]
+        scatStrength = stack.layers[0].scatterers.strength[0]
+
+        pOut = stack.backscatter(self.pIn)
+
+        self.assertAlmostEqual(pOut.Ex, scatStrength * exp(1j * self.pIn.k * 2 * scatPosition))
+        self.assertAlmostEqual(pOut.Ey, 0)
+
+    def testBackscatterSingleLayerNoBirefringence(self):
+        stack = TissueStackSingleLayerNoBirefringence()
+        expectedEx = 0
+
+        scatterers = stack.layers[0].scatterers
+        for (dz, strength) in zip(scatterers.dz, scatterers.strength):
+            scatPosition = stack.offset + dz
+            expectedEx += strength * exp(1j * self.pIn.k * 2 * scatPosition)
+
+        pOut = stack.backscatter(self.pIn)
+
+        self.assertAlmostEqual(pOut.Ex, expectedEx)
+        self.assertAlmostEqual(pOut.Ey, 0)
+
+    def testBackscatterSingleLayerPerpendicular(self):
+        """ With a beam polarization perpendicular to the one and only tissue optic axis, the retarding effect
+         of each scatterer is simply ikL(1+dn). """
+        self.pIn = JonesVector.vertical()
+        self.pIn.k = self.k
+        stack = TissueStackSingleHorizontal()
+        layer = stack.layers[0]
+
+        expectedEy = 0
+        for (dz, strength) in zip(layer.scatterers.dz, layer.scatterers.strength):
+            expectedEy += strength * exp(1j * self.pIn.k * 2 * (stack.offset + dz + dz * layer.birefringence))
 
-        # no ref to compare with
+        pOut = stack.backscatter(self.pIn)
+
+        self.assertAlmostEqual(pOut.Ex, 0)
+        self.assertAlmostEqual(pOut.Ey, expectedEy)
+
+    def testBackscatterSingleLayer(self):
+        stack = TissueStackSingle()
+        layer = stack.layers[0]
+
+        expectedEx = 0
+        expectedEy = 0
+        for (dz, strength) in zip(layer.scatterers.dz, layer.scatterers.strength):
+            jPhi = 1j * self.pIn.k * 2 * dz
+            A = exp(jPhi)
+            D = exp(jPhi * (1 + layer.birefringence))
+            a = layer.orientation
+            Ap = cos(a) ** 2 * A + D * sin(a) ** 2
+            Bp = (A - D) * sin(a) * cos(a)
+            Dp = cos(a) ** 2 * D + A * sin(a) ** 2
+
+            expectedEx += Ap * strength * exp(1j * self.pIn.k * 2 * stack.offset)
+            expectedEy += Bp * strength * exp(1j * self.pIn.k * 2 * stack.offset)
+
+        pOut = stack.backscatter(self.pIn)
+
+        self.assertAlmostEqual(pOut.Ex, expectedEx)
+        self.assertAlmostEqual(pOut.Ey, expectedEy)
+
+    def testBackscatter(self):
+        # combine testPropagateThrough with testBackscatterSingleLayer
+        # todo
+        pass
 
     def testBackscatterMany(self):
+        """ Qualitative test for pulse backscatter (multiple vectors). """
+        # todo: quantitative test
+
+        stack = TissueStackUnit()
         res = 5
         pIn = Pulse.horizontal(centerWavelength=1.3, wavelengthBandwidth=0.13, resolution=5)
 
-        pOut = self.stack.backscatterMany(pIn)
+        pOut = stack.backscatterMany(pIn)
 
         self.assertTrue(len(pOut) == res)
-        self.assertTrue(pOut[0].orientation != pOut[res//2].orientation)
+        self.assertTrue(pOut[0].orientation != pOut[res // 2].orientation)
 
 
-class TissueStackTestUnit(TissueStack):
+class TissueStackUnit(TissueStack):
     def __init__(self):
         layers = [TissueLayer(birefringence=0.001, opticAxis=(1, 0, 0), scattDensity=20, thickness=400),
                   TissueLayer(birefringence=0.002, opticAxis=(1, 1, 0), scattDensity=10, thickness=600),
                   TissueLayer(birefringence=0.003, opticAxis=(0, 1, 0), scattDensity=20, thickness=800)]
-        super(TissueStackTestUnit, self).__init__(offset=100, layers=layers)
+        super(TissueStackUnit, self).__init__(offset=100, layers=layers)
+
+
+class TissueStackNoBirefringence(TissueStack):
+    def __init__(self):
+        layers = [TissueLayer(birefringence=0, opticAxis=(1, 0, 0), scattDensity=20, thickness=400),
+                  TissueLayer(birefringence=0, opticAxis=(1, 1, 0), scattDensity=10, thickness=600),
+                  TissueLayer(birefringence=0, opticAxis=(0, 1, 0), scattDensity=20, thickness=800)]
+        super(TissueStackNoBirefringence, self).__init__(offset=100, layers=layers)
+
+
+class TissueStackOrientedHorizontally(TissueStack):
+    def __init__(self):
+        layers = [TissueLayer(birefringence=0.001, opticAxis=(1, 0, 0), scattDensity=20, thickness=400),
+                  TissueLayer(birefringence=0.002, opticAxis=(1, 0, 0), scattDensity=10, thickness=600),
+                  TissueLayer(birefringence=0.003, opticAxis=(1, 0, 0), scattDensity=20, thickness=800)]
+        super(TissueStackOrientedHorizontally, self).__init__(offset=100, layers=layers)
+
+
+class TissueStackOriented45Degrees(TissueStack):
+    def __init__(self):
+        layers = [TissueLayer(birefringence=0.001, opticAxis=(0, 1, 0), scattDensity=20, thickness=400),
+                  TissueLayer(birefringence=0.002, opticAxis=(0, 1, 0), scattDensity=10, thickness=600),
+                  TissueLayer(birefringence=0.003, opticAxis=(0, 1, 0), scattDensity=20, thickness=800)]
+        super(TissueStackOriented45Degrees, self).__init__(offset=100, layers=layers)
+
+
+class TissueStackSingleHorizontal(TissueStack):
+    def __init__(self):
+        layers = [TissueLayer(birefringence=0.001, opticAxis=(1, 0, 0), scattDensity=0.05, thickness=400)]
+        super(TissueStackSingleHorizontal, self).__init__(offset=100, layers=layers)
+
+
+class TissueStackSingle(TissueStack):
+    def __init__(self):
+        layers = [TissueLayer(birefringence=0.001, opticAxis=(0.6, 1.1, 0), scattDensity=20, thickness=400)]
+        super(TissueStackSingle, self).__init__(offset=100, layers=layers)
+
+
+class TissueStackSingleScattererSingleLayerNoBirefringence(TissueStack):
+    def __init__(self):
+        layers = [TissueLayer(birefringence=0, opticAxis=(1, 0, 0), scattDensity=20, thickness=400)]
+        for layer in layers:
+            layer.scatterers = Scatterers(layer.thickness, N=1)
+        super(TissueStackSingleScattererSingleLayerNoBirefringence, self).__init__(offset=100, layers=layers)
+
+
+class TissueStackSingleLayerNoBirefringence(TissueStack):
+    def __init__(self):
+        layers = [TissueLayer(birefringence=0, opticAxis=(1, 0, 0), scattDensity=0.05, thickness=400)]
+        super(TissueStackSingleLayerNoBirefringence, self).__init__(offset=100, layers=layers)
 
 
 if __name__ == '__main__':
diff --git a/polarization/tests/testsVector.py b/polarization/tests/testsVector.py
index 9a8253b..3186af9 100644
--- a/polarization/tests/testsVector.py
+++ b/polarization/tests/testsVector.py
@@ -1,32 +1,33 @@
-import envtest # modifies path
+import envtest  # modifies path
 from polarization.vector import *
 import numpy as np
 import random
 
 inf = float("+inf")
 
+
 class TestVector(envtest.MyTestCase):
     def randomVector(self):
-        x = random.random()*2-1
-        y = random.random()*2-1
-        z = random.random()*2-1
-        return Vector(x,y,z)
+        x = random.random() * 2 - 1
+        y = random.random() * 2 - 1
+        z = random.random() * 2 - 1
+        return Vector(x, y, z)
 
     def randomVectors(self, N):
         vectors = []
         for i in range(N):
-            x = random.random()*2-1
-            y = random.random()*2-1
-            z = random.random()*2-1
-            vectors.append( Vector(x,y,z) )
+            x = random.random() * 2 - 1
+            y = random.random() * 2 - 1
+            z = random.random() * 2 - 1
+            vectors.append(Vector(x, y, z))
 
     def randomUnitVectors(self, N):
         vectors = []
         for i in range(N):
-            x = random.random()*2-1
-            y = random.random()*2-1
-            z = random.random()*2-1
-            vectors.append( UnitVector(x,y,z) )
+            x = random.random() * 2 - 1
+            y = random.random() * 2 - 1
+            z = random.random() * 2 - 1
+            vectors.append(UnitVector(x, y, z))
 
     def testNullVector(self):
         v = Vector()
@@ -36,7 +37,7 @@ def testNullVector(self):
         self.assertEqual(v.z, 0)
 
     def testGetItemVector(self):
-        v = Vector(1,2,3)
+        v = Vector(1, 2, 3)
         self.assertIsNotNone(v)
         self.assertEqual(v[0], 1)
         self.assertEqual(v[1], 2)
@@ -45,14 +46,14 @@ def testGetItemVector(self):
             v[3]
 
     def testVectorValueInit(self):
-        v = Vector(1,2,3)
+        v = Vector(1, 2, 3)
         self.assertIsNotNone(v)
         self.assertEqual(v.x, 1)
         self.assertEqual(v.y, 2)
         self.assertEqual(v.z, 3)
 
     def testVectorCopy(self):
-        constVector = Vector(1,2,3)
+        constVector = Vector(1, 2, 3)
         v = Vector(constVector)
         self.assertIsNotNone(v)
         self.assertEqual(v.x, 1)
@@ -67,24 +68,24 @@ def testInvalidInput(self):
             Vector("ouch")
 
     def testConstVectors(self):
-        self.assertTrue(oHat.isEqualTo(Vector(0,0,0)))
-        self.assertTrue(xHat.isEqualTo(Vector(1,0,0)))
-        self.assertTrue(yHat.isEqualTo(Vector(0,1,0)))
-        self.assertTrue(zHat.isEqualTo(Vector(0,0,1)))
+        self.assertTrue(oHat.isEqualTo(Vector(0, 0, 0)))
+        self.assertTrue(xHat.isEqualTo(Vector(1, 0, 0)))
+        self.assertTrue(yHat.isEqualTo(Vector(0, 1, 0)))
+        self.assertTrue(zHat.isEqualTo(Vector(0, 0, 1)))
 
     def testEqualVectors(self):
-        self.assertEqual(oHat, Vector(0,0,0))
-        self.assertEqual(xHat, Vector(1,0,0))
-        self.assertEqual(yHat, Vector(0,1,0))
-        self.assertEqual(zHat, Vector(0,0,1))
+        self.assertEqual(oHat, Vector(0, 0, 0))
+        self.assertEqual(xHat, Vector(1, 0, 0))
+        self.assertEqual(yHat, Vector(0, 1, 0))
+        self.assertEqual(zHat, Vector(0, 0, 1))
 
-        self.assertNotEqual(Vector(1,2,3), Vector(0,2,3))
-        self.assertNotEqual(Vector(1,2,3), Vector(1,0,3))
-        self.assertNotEqual(Vector(1,2,3), Vector(1,2,0))
+        self.assertNotEqual(Vector(1, 2, 3), Vector(0, 2, 3))
+        self.assertNotEqual(Vector(1, 2, 3), Vector(1, 0, 3))
+        self.assertNotEqual(Vector(1, 2, 3), Vector(1, 2, 0))
 
-        self.assertFalse(Vector(1,2,3).isAlmostEqualTo(Vector(1.1,2,3)))
-        self.assertFalse(Vector(1,2,3).isAlmostEqualTo(Vector(1,2.1,3)))
-        self.assertFalse(Vector(1,2,3).isAlmostEqualTo(Vector(1,2,3.1)))
+        self.assertFalse(Vector(1, 2, 3).isAlmostEqualTo(Vector(1.1, 2, 3)))
+        self.assertFalse(Vector(1, 2, 3).isAlmostEqualTo(Vector(1, 2.1, 3)))
+        self.assertFalse(Vector(1, 2, 3).isAlmostEqualTo(Vector(1, 2, 3.1)))
 
     def testCanSetValuesVector(self):
         v = Vector()
@@ -95,8 +96,8 @@ def testCanSetValuesVector(self):
         self.assertEqual(v.y, 2)
         self.assertEqual(v.z, 3)
 
-    def testCannotSetConstVectors(self):      
-        v = ConstVector(1,2,3)
+    def testCannotSetConstVectors(self):
+        v = ConstVector(1, 2, 3)
 
         with self.assertRaises(RuntimeError):
             v.x = 1
@@ -126,8 +127,8 @@ def testCannotSetConstVectors(self):
         with self.assertRaises(RuntimeError):
             zHat.z = 2
 
-    def testCannotNormalizeConstVectors(self):      
-        v = ConstVector(1,2,3)
+    def testCannotNormalizeConstVectors(self):
+        v = ConstVector(1, 2, 3)
 
         with self.assertRaises(RuntimeError):
             v.normalize()
@@ -141,51 +142,51 @@ def testVectorIsUnitary(self):
         self.assertTrue(zHat.isUnitary)
 
     def testNormalizedCrossedProduct(self):
-        v1 = Vector(1,2,3)
+        v1 = Vector(1, 2, 3)
         v2 = v1.anyPerpendicular()
         self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), 1.0, 6)
 
-        v1 = Vector(1,0,0)
-        v2 = Vector(1,1,0)
-        self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(np.pi/4), 6)
+        v1 = Vector(1, 0, 0)
+        v2 = Vector(1, 1, 0)
+        self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(np.pi / 4), 6)
 
-        v1 = Vector(2,0,0)
-        v2 = Vector(2,2,0)
-        self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(np.pi/4), 6)
+        v1 = Vector(2, 0, 0)
+        v2 = Vector(2, 2, 0)
+        self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(np.pi / 4), 6)
 
-        v1 = Vector(0,0,1)
-        v2 = Vector(0,1,1)
-        self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(np.pi/4), 6)
+        v1 = Vector(0, 0, 1)
+        v2 = Vector(0, 1, 1)
+        self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(np.pi / 4), 6)
 
-        v1 = Vector(0,1,0)
-        v2 = Vector(0,1,1)
-        self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(np.pi/4), 6)
+        v1 = Vector(0, 1, 0)
+        v2 = Vector(0, 1, 1)
+        self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(np.pi / 4), 6)
 
-        v1 = Vector(0,1,0)
-        v2 = Vector(0,0,0)
+        v1 = Vector(0, 1, 0)
+        v2 = Vector(0, 0, 0)
         self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), 0, 6)
 
     def testOperators(self):
-        v1 = Vector(1,2,3)
-        self.assertTrue( (v1*2).isAlmostEqualTo(Vector(2,4,6)))
-        self.assertTrue( (2*v1).isAlmostEqualTo(Vector(2,4,6)))
-        self.assertTrue( (v1/2.0).isAlmostEqualTo(Vector(0.5,1,1.5)))
-        self.assertTrue( (v1+v1+v1).isAlmostEqualTo(Vector(3,6,9)))
-        self.assertTrue( (v1-2*v1).isAlmostEqualTo(Vector(-1,-2,-3)))
-        self.assertTrue( (-v1).isAlmostEqualTo(Vector(-1,-2,-3)))
+        v1 = Vector(1, 2, 3)
+        self.assertTrue((v1 * 2).isAlmostEqualTo(Vector(2, 4, 6)))
+        self.assertTrue((2 * v1).isAlmostEqualTo(Vector(2, 4, 6)))
+        self.assertTrue((v1 / 2.0).isAlmostEqualTo(Vector(0.5, 1, 1.5)))
+        self.assertTrue((v1 + v1 + v1).isAlmostEqualTo(Vector(3, 6, 9)))
+        self.assertTrue((v1 - 2 * v1).isAlmostEqualTo(Vector(-1, -2, -3)))
+        self.assertTrue((-v1).isAlmostEqualTo(Vector(-1, -2, -3)))
 
     def testRotateAround(self):
-        v1 = Vector(1,2,3)
+        v1 = Vector(1, 2, 3)
         axis = v1.anyUnitaryPerpendicular()
         v1.rotateAround(axis, 0.5)
-        v2 = Vector(1,2,3)
+        v2 = Vector(1, 2, 3)
         self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), np.sin(0.5), 6)
 
         for i in range(1000):
             v1 = self.randomVector()
             v2 = Vector(v1)
             axis = v1.anyUnitaryPerpendicular()
-            angle = 2*np.pi*random.random()
+            angle = 2 * np.pi * random.random()
             v1.rotateAround(axis, angle)
             self.assertAlmostEqual(v1.normalizedCrossProduct(v2).abs(), abs(np.sin(angle)), 6)
 
@@ -194,30 +195,30 @@ def testAngleWith(self):
             v1 = self.randomVector()
             v2 = Vector(v1)
             axis = v1.anyUnitaryPerpendicular()
-            angle = np.pi*(random.random()*2-1)
+            angle = np.pi * (random.random() * 2 - 1)
             v1.rotateAround(axis, angle)
             self.assertAlmostEqual(v1.angleWith(v2, axis), -angle, 6)
 
-    def testVectorIsParallel(self):     
-        v1 = Vector(1,2,3)
-        v2 = Vector(2,4,6)
+    def testVectorIsParallel(self):
+        v1 = Vector(1, 2, 3)
+        v2 = Vector(2, 4, 6)
 
         self.assertTrue(v1.isParallelTo(v2))
         self.assertFalse(v1.isParallelTo(xHat))
         self.assertFalse(v1.isParallelTo(yHat))
         self.assertFalse(v1.isParallelTo(zHat))
 
-    def testVectorIsAlmostParallel(self):     
-        v1 = Vector(1,2,3)
-        v3 = Vector(1,2,3.0000001)
+    def testVectorIsAlmostParallel(self):
+        v1 = Vector(1, 2, 3)
+        v3 = Vector(1, 2, 3.0000001)
         self.assertTrue(v1.isParallelTo(v3))
 
-    def testVectorIsAlmostNotParallel(self):     
-        v1 = Vector(1,2,3)
-        v3 = Vector(1,2,3.0001)
+    def testVectorIsAlmostNotParallel(self):
+        v1 = Vector(1, 2, 3)
+        v3 = Vector(1, 2, 3.0001)
         self.assertFalse(v1.isParallelTo(v3))
 
-    def testVectorIsPerpendicular(self):     
+    def testVectorIsPerpendicular(self):
         self.assertTrue(xHat.isPerpendicularTo(yHat))
         self.assertTrue(xHat.isPerpendicularTo(zHat))
         self.assertFalse(xHat.isPerpendicularTo(xHat))
@@ -232,15 +233,15 @@ def testVectorIsPerpendicular(self):
 
     def testManyPerpendicularVectors(self):
         for i in range(10000):
-            v1 = Vector(random.random()*2-1, random.random()*2-1, random.random()*2-1)
-            v2 = Vector(random.random()*2-1, random.random()*2-1, random.random()*2-1)
+            v1 = Vector(random.random() * 2 - 1, random.random() * 2 - 1, random.random() * 2 - 1)
+            v2 = Vector(random.random() * 2 - 1, random.random() * 2 - 1, random.random() * 2 - 1)
             self.assertTrue(v1.cross(v2).isPerpendicularTo(v1))
             self.assertTrue(v1.cross(v2).isPerpendicularTo(v2))
 
     def testAngle(self):
-        self.assertEqual(xHat.angleWith(yHat, zHat), np.pi/2)
-        self.assertEqual(yHat.angleWith(xHat, zHat), -np.pi/2)
-        self.assertEqual(xHat.angleWith(yHat, -zHat),-np.pi/2)
+        self.assertEqual(xHat.angleWith(yHat, zHat), np.pi / 2)
+        self.assertEqual(yHat.angleWith(xHat, zHat), -np.pi / 2)
+        self.assertEqual(xHat.angleWith(yHat, -zHat), -np.pi / 2)
 
     def testQuadrantAngle(self):
         plus45 = xHat + yHat
@@ -248,30 +249,30 @@ def testQuadrantAngle(self):
         minus45 = xHat - yHat
         minus135 = -xHat - yHat
 
-        self.assertAlmostEqual(xHat.angleWith(plus45, zHat), np.pi/4, 6)
-        self.assertAlmostEqual(xHat.angleWith(plus135, zHat), 3*np.pi/4, 6)
+        self.assertAlmostEqual(xHat.angleWith(plus45, zHat), np.pi / 4, 6)
+        self.assertAlmostEqual(xHat.angleWith(plus135, zHat), 3 * np.pi / 4, 6)
 
-        self.assertAlmostEqual(xHat.angleWith(minus45, zHat), -np.pi/4, 6)
-        self.assertAlmostEqual(xHat.angleWith(minus135, zHat), -3*np.pi/4, 6)
+        self.assertAlmostEqual(xHat.angleWith(minus45, zHat), -np.pi / 4, 6)
+        self.assertAlmostEqual(xHat.angleWith(minus135, zHat), -3 * np.pi / 4, 6)
 
-        self.assertAlmostEqual(xHat.angleWith(plus45, -zHat), -np.pi/4, 6)
-        self.assertAlmostEqual(xHat.angleWith(plus135, -zHat), -3*np.pi/4, 6)
+        self.assertAlmostEqual(xHat.angleWith(plus45, -zHat), -np.pi / 4, 6)
+        self.assertAlmostEqual(xHat.angleWith(plus135, -zHat), -3 * np.pi / 4, 6)
 
-        self.assertAlmostEqual(xHat.angleWith(minus45, -zHat), np.pi/4, 6)
-        self.assertAlmostEqual(xHat.angleWith(minus135, -zHat), 3*np.pi/4, 6)
+        self.assertAlmostEqual(xHat.angleWith(minus45, -zHat), np.pi / 4, 6)
+        self.assertAlmostEqual(xHat.angleWith(minus135, -zHat), 3 * np.pi / 4, 6)
 
         self.assertAlmostEqual(xHat.angleWith(-xHat, zHat), -np.pi, 6)
 
     def testScaledSum(self):
-        v1 = Vector(1,2,3)
-        v2 = Vector(4,5,6)
-        s  = 2
+        v1 = Vector(1, 2, 3)
+        v2 = Vector(4, 5, 6)
+        s = 2
         v3 = Vector(9, 12, 15)
         self.assertEqual(Vector.fromScaledSum(v1, v2, s), v3)
         self.assertEqual(v1.addScaled(v2, s), v3)
 
     def testAnyPerpendicular(self):
-        vectors = (Vector(1,2,3), Vector(-1,-2-3), xHat, yHat, zHat )
+        vectors = (Vector(1, 2, 3), Vector(-1, -2 - 3), xHat, yHat, zHat)
         for v in vectors:
             self.assertTrue(v.anyPerpendicular().isPerpendicularTo(v))
 
@@ -279,7 +280,7 @@ def testAnyPerpendicular(self):
         self.assertIsNone(oHat.anyPerpendicular())
 
     def testIsInKnownPlane(self):
-        v = Vector(1,2,3)
+        v = Vector(1, 2, 3)
         self.assertTrue(v.isInXYPlane(atZ=3))
         self.assertTrue(v.isInYZPlane(atX=1))
         self.assertTrue(v.isInZXPlane(atY=2))
@@ -289,57 +290,56 @@ def testIsInKnownPlane(self):
         self.assertFalse(v.isInZXPlane(atY=2.1))
 
     def testIsInSomePlane(self):
-        v = Vector(1,1,0)
-        self.assertTrue(v.isInPlane(origin=Vector(0,0,0), normal=zHat))
-        self.assertFalse(v.isInPlane(origin=Vector(1,1,1), normal=zHat))
+        v = Vector(1, 1, 0)
+        self.assertTrue(v.isInPlane(origin=Vector(0, 0, 0), normal=zHat))
+        self.assertFalse(v.isInPlane(origin=Vector(1, 1, 1), normal=zHat))
 
     def testRotationsUnitVectors(self):
         for i in range(10000):
-            v1 = UnitVector(random.random()*2-1, random.random()*2-1, random.random()*2-1)
-            v2 = UnitVector(random.random()*2-1, random.random()*2-1, random.random()*2-1)
+            v1 = UnitVector(random.random() * 2 - 1, random.random() * 2 - 1, random.random() * 2 - 1)
+            v2 = UnitVector(random.random() * 2 - 1, random.random() * 2 - 1, random.random() * 2 - 1)
             axis = v1.cross(v2)
             angle = v1.angleWith(v2, axis)
             v1.rotateAround(axis, angle)
 
-            self.assertTrue( (v1-v2).abs() < 1e-6)
-            self.assertTrue( v1.isUnitary)
-            self.assertTrue( v2.isUnitary)
+            self.assertTrue((v1 - v2).abs() < 1e-6)
+            self.assertTrue(v1.isUnitary)
+            self.assertTrue(v2.isUnitary)
 
-        self.assertTrue(Vector(1,0,0).rotateAround(zHat, np.pi/2).isAlmostEqualTo( yHat))
-        self.assertTrue(Vector(0,1,0).rotateAround(xHat, np.pi/2).isAlmostEqualTo( zHat))
-        self.assertTrue(Vector(0,0,1).rotateAround(yHat, np.pi/2).isAlmostEqualTo( xHat))
-            
+        self.assertTrue(Vector(1, 0, 0).rotateAround(zHat, np.pi / 2).isAlmostEqualTo(yHat))
+        self.assertTrue(Vector(0, 1, 0).rotateAround(xHat, np.pi / 2).isAlmostEqualTo(zHat))
+        self.assertTrue(Vector(0, 0, 1).rotateAround(yHat, np.pi / 2).isAlmostEqualTo(xHat))
 
     def testRotationsVectors(self):
         for i in range(10000):
-            v1 = Vector(random.random()*2-1, random.random()*2-1, random.random()*2-1)
-            v2 = Vector(random.random()*2-1, random.random()*2-1, random.random()*2-1)
+            v1 = Vector(random.random() * 2 - 1, random.random() * 2 - 1, random.random() * 2 - 1)
+            v2 = Vector(random.random() * 2 - 1, random.random() * 2 - 1, random.random() * 2 - 1)
             axis = v1.cross(v2)
             angle = v1.angleWith(v2, axis)
             v1.rotateAround(axis, angle)
 
-            self.assertTrue( v1.dot(v2) >= 0)
-            self.assertTrue( v1.cross(v2).abs() < 1e-6)
-            self.assertTrue( v1.isParallelTo(v2) )
+            self.assertTrue(v1.dot(v2) >= 0)
+            self.assertTrue(v1.cross(v2).abs() < 1e-6)
+            self.assertTrue(v1.isParallelTo(v2))
 
     def testPlaneOfIncidence(self):
-        ez = Vector(0,1,1).normalized()
+        ez = Vector(0, 1, 1).normalized()
         plane = ez.planeOfIncidence(normal=zHat)
         self.assertTrue(plane.isParallelTo(xHat))
         self.assertTrue(plane.isUnitary)
 
-        ez = Vector(0,1,1).normalized()
+        ez = Vector(0, 1, 1).normalized()
         plane = ez.planeOfIncidence(normal=-zHat)
         self.assertTrue(plane.isParallelTo(xHat))
         self.assertTrue(plane.isUnitary)
 
-        ez = Vector(0,-1,1).normalized()
+        ez = Vector(0, -1, 1).normalized()
         plane = ez.planeOfIncidence(normal=zHat)
         self.assertTrue(plane.isParallelTo(xHat))
         self.assertTrue(plane.isUnitary)
 
-        ez = Vector(0,1,1) # unnormalized
-        plane = ez.planeOfIncidence(normal=Vector(0,0,10))
+        ez = Vector(0, 1, 1)  # unnormalized
+        plane = ez.planeOfIncidence(normal=Vector(0, 0, 10))
         self.assertTrue(plane.isParallelTo(xHat))
         self.assertTrue(plane.isUnitary)
 
@@ -356,42 +356,42 @@ def testPlaneOfIncidence(self):
         self.assertTrue(plane.isUnitary)
 
     def testAngleOfIncidence(self):
-        ez = Vector(0,1,1).normalized()
-        surfaceNormal = Vector(0,0,1)
+        ez = Vector(0, 1, 1).normalized()
+        surfaceNormal = Vector(0, 0, 1)
         angle, planeNormal, actualNormal = ez.angleOfIncidence(normal=surfaceNormal)
-        self.assertAlmostEqual(ez.angleWith(surfaceNormal, planeNormal), np.pi/4,6)
+        self.assertAlmostEqual(ez.angleWith(surfaceNormal, planeNormal), np.pi / 4, 6)
 
-        surfaceNormal = Vector(0,0,-1)
+        surfaceNormal = Vector(0, 0, -1)
         angle, planeNormal, actualNormal = ez.angleOfIncidence(normal=surfaceNormal)
-        self.assertAlmostEqual(ez.angleWith(actualNormal, planeNormal), np.pi/4,6)
+        self.assertAlmostEqual(ez.angleWith(actualNormal, planeNormal), np.pi / 4, 6)
 
-        ez = Vector(0,-1,1).normalized()
+        ez = Vector(0, -1, 1).normalized()
         surfaceNormal = zHat
         angle, planeNormal, actualNormal = ez.angleOfIncidence(normal=surfaceNormal)
-        self.assertAlmostEqual(ez.angleWith(actualNormal, planeNormal), np.pi/4,6)
+        self.assertAlmostEqual(ez.angleWith(actualNormal, planeNormal), np.pi / 4, 6)
 
-        ez = Vector(0,1,1).normalized()
+        ez = Vector(0, 1, 1).normalized()
         surfaceNormal = zHat
         angle, planeNormal, actualNormal = ez.angleOfIncidence(surfaceNormal)
-        self.assertAlmostEqual(angle, np.pi/4,6)
+        self.assertAlmostEqual(angle, np.pi / 4, 6)
         self.assertTrue(planeNormal.isUnitary)
 
-        ez = Vector(0,-1,1).normalized()
+        ez = Vector(0, -1, 1).normalized()
         surfaceNormal = zHat
         angle, planeNormal, actualNormal = ez.angleOfIncidence(surfaceNormal)
-        self.assertAlmostEqual(angle, np.pi/4,6)
+        self.assertAlmostEqual(angle, np.pi / 4, 6)
         self.assertTrue(planeNormal.isUnitary)
 
         ez = zHat
         surfaceNormal = zHat
         angle, planeNormal, actualNormal = ez.angleOfIncidence(surfaceNormal)
-        self.assertAlmostEqual(angle, 0,6)
+        self.assertAlmostEqual(angle, 0, 6)
         self.assertTrue(planeNormal.isUnitary)
 
-        ez = Vector(0,0,1)
-        surfaceNormal = Vector(0,0,1)
+        ez = Vector(0, 0, 1)
+        surfaceNormal = Vector(0, 0, 1)
         angle, planeNormal, actualNormal = ez.angleOfIncidence(surfaceNormal)
-        self.assertAlmostEqual(angle, 0,6)
+        self.assertAlmostEqual(angle, 0, 6)
         self.assertTrue(planeNormal.isUnitary)