diff --git a/examples/test_LRM.py b/examples/dyakonov_nikitin_2009.py
similarity index 59%
rename from examples/test_LRM.py
rename to examples/dyakonov_nikitin_2009.py
index ff15418..fb5cc67 100644
--- a/examples/test_LRM.py
+++ b/examples/dyakonov_nikitin_2009.py
@@ -1,16 +1,18 @@
-__author__ = 'Pavel Dmitriev'
-
 import numpy as np
 import matplotlib.pyplot as plt
-from uniaxialTransferMatrix import *
-from isotropicTransferMatrix import *
-from transferMatrix import solve_transfer_matrix
-import scipy.linalg
+
+from PyATMM.uniaxialTransferMatrix import *
+from PyATMM.isotropicTransferMatrix import *
+from PyATMM.transferMatrix import solve_transfer_matrix
+
+__author__ = 'Pavel Dmitriev'
+
+#
+# Setup
+#
 
 c = 299792458.  # m/c
 w = c
-#w = 2*np.pi * 500. * 10 **12
-
 
 #
 # Theta sweep
@@ -29,7 +31,6 @@
 #
 eps_Ta2O5 = 4.41
 ran_d = np.linspace(5, 0, 50, endpoint=True)
-#ran_d = [3.]
 
 #
 # Anisotropic layer, YV04
@@ -45,8 +46,6 @@
 
 a_refl_ps = []
 a_refl_sp = []
-
-
 for d in ran_d:
     a_refl_pp_ = []
     a_refl_ss_ = []
@@ -55,18 +54,22 @@
     a_refl_sp_ = []
 
     d *= 2.*np.pi
-    #d *= 2 * np.pi * (c / w)
 
     for kx in ran_kx:
+        # Layer 0
         D_ZnSe = build_isotropic_dynamic_matrix(eps_ZnSe, w, kx, ky)
 
+        # Layer 1
         D_Ta2O5 = build_isotropic_dynamic_matrix(eps_Ta2O5, w, kx, ky)
         D_Ta2O5_p = build_isotropic_propagation_matrix(eps_Ta2O5, w, kx, ky, d)
+
+        # Layer 2
         D_YVo4 = build_uniaxial_dynamic_matrix(eps_YVO4_per, eps_YVO4_par, w, kx, ky, opticAxis=[np.cos(phi), np.sin(phi), 0])
 
+        # Multiply Matricies
         D = np.dot(np.dot(np.dot(np.linalg.inv(D_ZnSe), D_Ta2O5), D_Ta2O5_p), np.dot(np.linalg.inv(D_Ta2O5), D_YVo4))
-        #D = np.dot(np.dot(np.linalg.inv(D_ZnSe), D_Ta2O5), D_Ta2O5_p)
 
+        # Solve Matrix
         r_ss, r_sp, r_ps, r_pp, t_ss, t_sp, t_ps, t_pp = solve_transfer_matrix(D)
         a_refl_pp_.append(np.abs(r_pp**2))
         a_refl_ss_.append(np.abs(r_ss**2))
@@ -80,40 +83,17 @@
     a_refl_sp.append(a_refl_sp_)
     a_refl_ps.append(a_refl_ps_)
 
+
 plt.figure(figsize=(9, 6))
-plt.imshow(a_refl_ps, interpolation='none', extent=[np.min(theta), np.max(theta), np.min(ran_d), np.max(ran_d)], aspect="auto",  cmap=plt.get_cmap("jet"))
+plt.imshow(a_refl_ps, interpolation='none', extent=[np.min(theta), np.max(theta), np.min(ran_d), np.max(ran_d)],
+           aspect="auto",  cmap=plt.get_cmap("jet"))
+
 plt.colorbar()
 plt.xlabel("$\Theta$, degrees")
-#plt.xticks([54.18, 54.19, 54.20, 54.21, 54.22])
 plt.xticks([54.20, 54.205, 54.21, 54.215, 54.22])
 plt.ticklabel_format(style = 'sci', useOffset=False)
 plt.ylabel("$\Delta / \lambda$")
 plt.title("R_ps")
 
-#plt.suptitle("Otto Geometry test", fontsize=40)
-plt.savefig("otto_test.pdf", bbox_inches='tight')
-plt.show()
+plt.show(block=True)
 
-
-
-# #TE
-# plt.plot(theta, a_refl_ss[0])
-# plt.plot(theta, a_refl_sp[0])
-# #plt.plot(theta, numpy.add(a_refl_ss, a_refl_sp))
-# #plt.legend(['SS', 'SP', 'Sum'], loc='best')
-#
-# #TM
-# plt.plot(theta, a_refl_pp[0])
-# plt.plot(theta, a_refl_ps[0])
-# #plt.plot(theta, numpy.add(a_refl_pp, a_refl_ps))
-# #plt.legend(['PP', 'PS', 'Sum'], loc='best')
-#
-# #plt.legend(['SS', 'SP', 'Sum', 'PP', 'PS', 'Sum'], loc='best')
-# plt.legend(['SS', 'SP', 'PP', 'PS'], loc='best')
-# plt.xlabel("$\Theta$, degrees")
-# #plt.xticks([54.18, 54.19, 54.20, 54.21, 54.22])
-# plt.ticklabel_format(style = 'sci', useOffset=False)
-# plt.ylabel("Reflectance")
-# plt.title("$\Delta / \lambda = 3$")
-# plt.savefig("otto_slice_noanis.pdf", bbox_inches='tight')
-# plt.show(block=True)
\ No newline at end of file
diff --git a/examples/test_Dyakonov_Takayama2011.py b/examples/dyakonov_takayama_2011.py
similarity index 100%
rename from examples/test_Dyakonov_Takayama2011.py
rename to examples/dyakonov_takayama_2011.py
diff --git a/examples/test_isotropic_2.py b/examples/isotropic_brewster.py
similarity index 70%
rename from examples/test_isotropic_2.py
rename to examples/isotropic_brewster.py
index 24e746e..e7990a3 100644
--- a/examples/test_isotropic_2.py
+++ b/examples/isotropic_brewster.py
@@ -1,12 +1,17 @@
-__author__ = 'Pavel Dmitriev'
-
 import numpy as np
 import matplotlib.pyplot as plt
-from isotropicTransferMatrix import *
-from transferMatrix import solve_transfer_matrix
-import PyTMM.pytmm.transferMatrix as tm
 import scipy.linalg
 
+import PyTMM.transferMatrix as tm
+
+from PyATMM.isotropicTransferMatrix import *
+from PyATMM.transferMatrix import solve_transfer_matrix
+
+__author__ = 'Pavel Dmitriev'
+
+#
+# Setup
+#
 c = 299792458  # m/c
 
 w = 500 * 10 **12
@@ -19,27 +24,28 @@
 
 ky = 0
 
+
+#
+# PyTMM
+#
 refl_p = []
 refl_s = []
 for kx in ran_kx:
     B = tm.TransferMatrix.boundingLayer(1, n_1, np.arcsin(kx*c/w), tm.Polarization.s)
     C = tm.TransferMatrix.boundingLayer(1, n_1, np.arcsin(kx*c/w), tm.Polarization.p)
-    #a.appendRight(b)
 
     M = scipy.linalg.block_diag(B.matrix, C.matrix)
     r_ss, r_sp, r_ps, r_pp, t_ss, t_sp, t_ps, t_pp = solve_transfer_matrix(M)
     refl_p.append(np.abs(r_pp**2))
     refl_s.append(np.abs(r_ss**2))
 
+#
+# PyATMM
+#
 a_refl_p = []
 a_refl_s = []
 for kx in ran_kx:
     D = build_isotropic_bounding_layer_matrix(eps_1, w, kx, ky)
-    #D = build_isotropic_layer_matrix(eps_1, w, kx, ky, d_m)
-    #D = build_isotropic_propagation_matrix(eps_1, w, kx, ky, d_m)
-
-    #r_ss, r_sp, r_ps, r_pp, t_ss, t_sp, t_ps, t_pp = solve_transfer_matrix(numpy.dot(D_b, D))
-    #a_refl.append(np.abs(r_pp**2))
 
     r_ss, r_sp, r_ps, r_pp, t_ss, t_sp, t_ps, t_pp = solve_transfer_matrix(D)
     a_refl_p.append(np.abs(r_pp**2))
@@ -50,16 +56,14 @@
 #PyATMM
 plt.plot(angle, a_refl_p)
 plt.plot(angle, a_refl_s)
-#plt.plot(ran_l*10**9, sum)
 
 #PyTMM
-#plt.plot(ran_l*10**9, refl0, 'o')
 plt.plot(angle, refl_p, '+')
 plt.plot(angle, refl_s, '+')
 
 
 plt.xlabel("kx, $m^{-1}$")
 plt.ylabel("Reflectance")
-#plt.title("Reflectance of ideal single-layer antireflective coating")
+plt.title("Brewster angle")
 plt.legend(['PP', 'SS', 'P', 'S'], loc='best')
 plt.show(block=True)
\ No newline at end of file
diff --git a/examples/test_isotropic.py b/examples/isotropic_slab.py
similarity index 83%
rename from examples/test_isotropic.py
rename to examples/isotropic_slab.py
index b6fd36d..a0ed5c6 100644
--- a/examples/test_isotropic.py
+++ b/examples/isotropic_slab.py
@@ -1,16 +1,23 @@
-__author__ = 'Pavel Dmitriev'
-
 import numpy as np
 import matplotlib.pyplot as plt
-from isotropicTransferMatrix import *
-from transferMatrix import solve_transfer_matrix
-import PyTMM.pytmm.transferMatrix as tm
 import scipy.linalg
 
+import PyTMM.transferMatrix as tm
+
+from PyATMM.isotropicTransferMatrix import *
+from PyATMM.transferMatrix import solve_transfer_matrix
+
+
+__author__ = 'Pavel Dmitriev'
+
+#
+# Setup
+#
+
 c = 299792458  # m/c
 
 ran_w = np.linspace(50, 1200, 100, endpoint=False)
-ran_w = np.multiply(ran_w, 10 ** 12)  # (1/(2*np.pi))*
+ran_w = np.multiply(ran_w, 10**12)  # (1/(2*np.pi))*
 ran_l = np.divide(c*(2*np.pi), ran_w)
 
 eps_1 = 2.25
@@ -22,6 +29,9 @@
 d_nm = 2000
 d_m = d_nm * 10**-9
 
+#
+# PyTMM
+#
 refl0 = []
 refl_p = []
 refl_s = []
@@ -34,10 +44,12 @@
     refl_p.append(np.abs(r_pp**2))
     refl_s.append(np.abs(r_ss**2))
 
+#
+# PyATMM
+#
 a_refl_p = []
 a_refl_s = []
 for w in ran_w:
-    #D = build_isotropic_bounding_layer_matrix(eps_1, w, kx, ky)
     D = build_isotropic_layer_matrix(eps_1, w, kx, ky, d_m)
 
     r_ss, r_sp, r_ps, r_pp, t_ss, t_sp, t_ps, t_pp = solve_transfer_matrix(D)
diff --git a/examples/test_isotropic_3.py b/examples/isotropic_totalinternal.py
similarity index 86%
rename from examples/test_isotropic_3.py
rename to examples/isotropic_totalinternal.py
index 2132eeb..bfd867c 100644
--- a/examples/test_isotropic_3.py
+++ b/examples/isotropic_totalinternal.py
@@ -1,12 +1,17 @@
-__author__ = 'Pavel Dmitriev'
-
 import numpy as np
 import matplotlib.pyplot as plt
-from isotropicTransferMatrix import *
-from transferMatrix import solve_transfer_matrix
-import PyTMM.pytmm.transferMatrix as tm
 import scipy.linalg
 
+import PyTMM.transferMatrix as tm
+
+from PyATMM.isotropicTransferMatrix import *
+from PyATMM.transferMatrix import solve_transfer_matrix
+
+__author__ = 'Pavel Dmitriev'
+
+#
+# Setup
+#
 c = 299792458  # m/c
 
 w = 500 * 10 **12
@@ -22,9 +27,12 @@
 
 ky = 0
 
+
+#
+# PyTMM
+#
 refl_p = []
 refl_s = []
-
 for kx in ran_kx:
     B = tm.TransferMatrix.boundingLayer(n_1, n_2, np.arcsin(kx*c/w), tm.Polarization.s)
     C = tm.TransferMatrix.boundingLayer(n_1, n_2, np.arcsin(kx*c/w), tm.Polarization.p)
@@ -35,6 +43,10 @@
     refl_p.append(np.abs(r_pp**2))
     refl_s.append(np.abs(r_ss**2))
 
+
+#
+# PyATMM
+#
 a_refl_p = []
 a_refl_s = []
 for kx in ran_kx:
@@ -55,6 +67,6 @@
 
 plt.xlabel("angle, degrees")
 plt.ylabel("Reflectance")
-#plt.title("Reflectance of ideal single-layer antireflective coating")
+plt.title("Total internal reflection")
 plt.legend(['PP', 'SS', 'P', 'S'], loc='best')
 plt.show(block=True)
\ No newline at end of file
diff --git a/examples/test.py b/examples/test.py
deleted file mode 100644
index 138d2e1..0000000
--- a/examples/test.py
+++ /dev/null
@@ -1,76 +0,0 @@
-import numpy as np
-from transferMatrix import *
-import scipy
-from scipy import linalg, matrix
-
-def null(A, eps=1e-15):
-    u, s, vh = scipy.linalg.svd(A)
-    null_mask = (s <= eps)
-    null_space = scipy.compress(null_mask, vh, axis=0)
-    return scipy.transpose(null_space)
-
-def build_matrix(kx, ky, kz, Eps, w):
-    c = 1 # m/c
-    return (w**2 / c**2 ) * Eps - numpy.dot([kx, ky, kz], [kx, ky, kz])*numpy.eye(3) + numpy.outer([kx, ky, kz], [kx, ky, kz])
-
-
-kx=0.01
-ky=0.02
-
-e1=1.02
-e2=1.01
-e3=1
-
-# theta=np.pi*3/5
-# phi=np.pi/4
-# psi=np.pi*7/4
-theta=0
-phi=0
-psi=0
-
-w=1
-d=5
-
-P_an, gamma_an, Eps = build_anisotropic_layer_matrix(e1=e1, e2=e2, e3=e3, theta=theta, phi=phi, psi=psi, w=w, kx=kx, ky=ky, d=d)
-
-# print('gamma 0', gamma_an[0])
-A = [build_matrix(kx=kx, ky=ky, kz=g, Eps=Eps, w=1) for g in gamma_an]
-# print(np.linalg.det(A))
-
-for a in A:
-    a[2, 0] = 0
-    a[2, 1] = 0
-    a[2, 2] = 0
-
-# print('A', A)
-#
-P_num = [null(a).T[0] for a in A]
-#
-# print("An", P_an[0])
-# print("Num", P_num)
-#
-#
-# #print(numpy.dot(numpy.dot(Eps, P_num), [kx, ky, gamma_an[0]]))
-#
-# print(numpy.dot(numpy.dot(Eps, P_an[0]), [kx, ky, gamma_an[0]]))
-
-# print(np.dot(A, P_num))
-#
-# A = build_matrix(kx=kx, ky=ky, kz=gamma_an[0], Eps=Eps, w=1)
-# print(np.dot(A, P_an[0]))
-# A = build_matrix(kx=kx, ky=ky, kz=gamma_an[1], Eps=Eps, w=1)
-# print(np.dot(A, P_an[1]))
-# A = build_matrix(kx=kx, ky=ky, kz=gamma_an[2], Eps=Eps, w=1)
-# print(np.dot(A, P_an[2]))
-# A = build_matrix(kx=kx, ky=ky, kz=gamma_an[3], Eps=Eps, w=1)
-# print(np.dot(A, P_an[3]))
-print("Null-space")
-print("Pol: ", np.abs(P_num[0]), "Kz: ", np.abs(gamma_an[0]))
-print("Pol: ", np.abs(P_num[1]), "Kz: ", np.abs(gamma_an[1]))
-print("Pol: ", np.abs(P_num[2]), "Kz: ", np.abs(gamma_an[2]))
-print("Pol: ", np.abs(P_num[3]), "Kz: ", np.abs(gamma_an[3]))
-print("Analytic")
-print("Pol: ", np.abs(P_an[0]), "Kz: ", np.abs(gamma_an[0]))
-print("Pol: ", np.abs(P_an[1]), "Kz: ", np.abs(gamma_an[1]))
-print("Pol: ", np.abs(P_an[2]), "Kz: ", np.abs(gamma_an[2]))
-print("Pol: ", np.abs(P_an[3]), "Kz: ", np.abs(gamma_an[3]))
\ No newline at end of file
diff --git a/examples/test_solve_quartic.py b/examples/test_solve_quartic.py
deleted file mode 100644
index 113ddde..0000000
--- a/examples/test_solve_quartic.py
+++ /dev/null
@@ -1,14 +0,0 @@
-from unittest import TestCase
-from util import *
-__author__ = 'p.dmitriev'
-
-
-class TestSolve_quartic(TestCase):
-    def test_solve_quartic(self):
-        a = 1.
-        b = 0.
-        c = 0.
-        d = 0.
-        e = -1.
-        print(solve_quartic(a, b, c, d, e))
-        self.fail()
diff --git a/examples/test_uniaxial_2.py b/examples/uniaxial_bukharev_2015_bsc_thesis.py
similarity index 78%
rename from examples/test_uniaxial_2.py
rename to examples/uniaxial_bukharev_2015_bsc_thesis.py
index a90d611..b08134a 100644
--- a/examples/test_uniaxial_2.py
+++ b/examples/uniaxial_bukharev_2015_bsc_thesis.py
@@ -1,17 +1,18 @@
-__author__ = 'Pavel Dmitriev'
-
 import numpy as np
 import matplotlib.pyplot as plt
-from uniaxialTransferMatrix import *
-from isotropicTransferMatrix import *
-from transferMatrix import solve_transfer_matrix
-import scipy.linalg
+from PyATMM.uniaxialTransferMatrix import *
+from PyATMM.isotropicTransferMatrix import *
+from PyATMM.transferMatrix import solve_transfer_matrix
+
+__author__ = 'Pavel Dmitriev'
 
+#
+# Setup
+#
 c = 299792458  # m/c
 
 w = 500 * 10 **12
 l = 2*np.pi*c / w
-#ran_kx = np.linspace(0, (w/c)*0.99, 100, endpoint=False)
 
 angle = np.linspace(0, 90, 500, endpoint=False)
 ran_kx = (w/c) * np.sin(np.deg2rad(angle))
@@ -37,7 +38,6 @@
 
 for kx in ran_kx:
     D_0 = build_isotropic_dynamic_matrix(eps_g, w, n_g*kx, ky)
-    #D_1 = build_isotropic_dynamic_matrix(eps_ex, w, n_g*kx, ky)
     D_1 = build_uniaxial_dynamic_matrix(eps_ord, eps_ex, w, n_g*kx, ky, opticAxis=([0, 1./np.sqrt(2), 1./np.sqrt(2)]))
     D = np.dot(np.linalg.inv(D_0), D_1)
 
@@ -64,5 +64,5 @@
 
 plt.xlabel("angle, degrees")
 plt.ylabel("Reflectance")
-#plt.title("Reflectance of ideal single-layer antireflective coating")
+plt.title("Uniaxial total internal reflection")
 plt.show(block=True)
\ No newline at end of file
diff --git a/examples/test_uniaxial.py b/examples/uniaxial_slab.py
similarity index 70%
rename from examples/test_uniaxial.py
rename to examples/uniaxial_slab.py
index db9167c..3de7b41 100644
--- a/examples/test_uniaxial.py
+++ b/examples/uniaxial_slab.py
@@ -1,12 +1,17 @@
-__author__ = 'Pavel Dmitriev'
-
 import numpy as np
 import matplotlib.pyplot as plt
-from uniaxialTransferMatrix import *
-from transferMatrix import solve_transfer_matrix
-import PyTMM.pytmm.transferMatrix as tm
 import scipy.linalg
 
+import PyTMM.transferMatrix as tm
+
+from PyATMM.uniaxialTransferMatrix import *
+from PyATMM.transferMatrix import solve_transfer_matrix
+
+__author__ = 'Pavel Dmitriev'
+
+#
+# Setup
+#
 c = 299792458  # m/c
 
 ran_w = np.linspace(50, 1200, 100, endpoint=False)
@@ -21,12 +26,12 @@
 kx = 0
 ky = 0
 
-# d_nm = 800 / (n_2 * 4) # quarter-wavelength coating
 d_nm = 2000
 d_m = d_nm * 10**-9
 
-
-
+#
+# PyTMM
+#
 refl_p = []
 refl_s = []
 for i in ran_l:
@@ -37,16 +42,13 @@
     refl_p.append(np.abs(r_pp**2))
     refl_s.append(np.abs(r_ss**2))
 
-
+#
+# PyATMM
+#
 a_refl_p = []
 a_refl_s = []
 for w in ran_w:
-    #D = build_isotropic_bounding_layer_matrix(eps_1, w, kx, ky)
     D = build_uniaxial_layer_matrix(eps_1, eps_2, w, kx, ky, d_m, opticAxis=([0., 1., 0.]))
-    #D = build_isotropic_propagation_matrix(eps_1, w, kx, ky, d_m)
-
-    #r_ss, r_sp, r_ps, r_pp, t_ss, t_sp, t_ps, t_pp = solve_transfer_matrix(numpy.dot(D_b, D))
-    #a_refl.append(np.abs(r_pp**2))
 
     r_ss, r_sp, r_ps, r_pp, t_ss, t_sp, t_ps, t_pp = solve_transfer_matrix(D)
     a_refl_p.append(np.abs(r_pp**2))
@@ -56,16 +58,14 @@
 #PyATMM
 plt.plot(ran_l*10**9, a_refl_p)
 plt.plot(ran_l*10**9, a_refl_s)
-#plt.plot(ran_l*10**9, sum)
 
 #PyTMM
 plt.plot(ran_l*10**9, refl_p, 'o')
 plt.plot(ran_l*10**9, refl_s, 'o')
 
 
-
 plt.xlabel("Wavelength, nm")
 plt.ylabel("Reflectance")
-plt.title("Reflectance of ideal single-layer antireflective coating")
+plt.title("Reflectance uniaxial crystal")
 plt.legend(['PP', 'SS'], loc='best')
 plt.show(block=True)
\ No newline at end of file