diff --git a/.github/workflows/pytest_and_coverage.yml b/.github/workflows/pytest_and_coverage.yml
index 6b061c8..b7a89a0 100644
--- a/.github/workflows/pytest_and_coverage.yml
+++ b/.github/workflows/pytest_and_coverage.yml
@@ -36,9 +36,9 @@ jobs:
       run: |
         python3 -m pip install flake8
         # stop the build if there are Python syntax errors or undefined names
-        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        flake8 ./composites --count --select=E9,F63,F7,F82 --show-source --statistics
         # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
-        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+        flake8 ./composites --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest and coverage report
       if: env.USE_COVERAGE == 'true'
       run: |
diff --git a/composites/__init__.py b/composites/__init__.py
index 0b1f8ea..c12b009 100644
--- a/composites/__init__.py
+++ b/composites/__init__.py
@@ -1,12 +1,12 @@
 r"""
-===================================================================
-Methods to calculate composite plate properties (:mod:`composites`)
-===================================================================
+==================================
+composites API (:mod:`composites`)
+==================================
 
 .. currentmodule::composites
 
-The ``composites`` module includes functions used to calculate plate properties
-for laminated composites and isotropic plates.
+The ``composites`` module includes functions used to calculate properties and
+perform analysis on laminated composites and isotropic plates.
 
 Classical and first-order shear deformation theories are supported. For
 classical plate theories or classical laminated plate theories (CLPT), the
@@ -54,6 +54,9 @@
 .. automodule:: composites.utils
     :members:
 
+.. automodule:: composites.kassapoglou
+    :members:
+
 """
 import os
 
diff --git a/composites/core.pxd b/composites/core.pxd
index 819031a..3e2475e 100644
--- a/composites/core.pxd
+++ b/composites/core.pxd
@@ -53,9 +53,10 @@ cdef class Laminate:
     cpdef void calc_scf(Laminate)
     cpdef void calc_equivalent_properties(Laminate)
     cpdef void calc_constitutive_matrix(Laminate)
-    cpdef void force_balanced(Laminate)
-    cpdef void force_orthotropic(Laminate)
-    cpdef void force_symmetric(Laminate)
+    cpdef void make_balanced(Laminate)
+    cpdef void make_orthotropic(Laminate)
+    cpdef void make_symmetric(Laminate)
+    cpdef void make_smeared(Laminate)
     cpdef LaminationParameters calc_lamination_parameters(Laminate)
 
 cdef class GradABDE:
diff --git a/composites/core.pyx b/composites/core.pyx
index 2124cf0..7fd216c 100644
--- a/composites/core.pyx
+++ b/composites/core.pyx
@@ -616,29 +616,31 @@ cdef class Laminate:
             self.E45 += ply.q45L*(hk - hk_1)
             self.E55 += ply.q55L*(hk - hk_1)
 
-    cpdef void force_balanced(Laminate self):
-        r"""Force a balanced laminate
+
+    cpdef void make_balanced(Laminate self):
+        r"""Make a balanced laminate
 
         The attributes `A_{16}`, `A_{26}`, `B_{16}`, `B_{26}` are set to zero
-        to force a balanced laminate.
+        to make a balanced laminate.
 
         """
         if self.offset != 0.:
-            raise RuntimeError('Laminates with offset cannot be forced balanced!')
+            raise RuntimeError('Laminates with offset cannot be made balanced!')
         self.A16 = 0.
         self.A26 = 0.
         self.B16 = 0.
         self.B26 = 0.
 
-    cpdef void force_orthotropic(Laminate self):
-        r"""Force an orthotropic laminate
+
+    cpdef void make_orthotropic(Laminate self):
+        r"""Make an orthotropic laminate
 
         The attributes `A_{16}`, `A_{26}`, `B_{16}`, `B_{26}`, `D_{16}`,
-        `D_{26}` are set to zero to force an orthotropic laminate.
+        `D_{26}` are set to zero to make an orthotropic laminate.
 
         """
         if self.offset != 0.:
-            raise RuntimeError('Laminates with offset cannot be forced orthotropic!')
+            raise RuntimeError('Laminates with offset cannot be made orthotropic!')
         self.A16 = 0.
         self.A26 = 0.
         self.B16 = 0.
@@ -646,15 +648,38 @@ cdef class Laminate:
         self.D16 = 0.
         self.D26 = 0.
 
-    cpdef void force_symmetric(Laminate self):
-        """Force a symmetric laminate
+
+    cpdef void make_symmetric(Laminate self):
+        """Make a symmetric laminate
 
         The `B_{ij}` terms of the constitutive matrix are set to zero.
 
         """
         if self.offset != 0.:
             raise RuntimeError(
-                    'Laminates with offset cannot be forced symmetric!')
+                    'Laminates with offset cannot be made symmetric!')
+        self.B11 = 0
+        self.B12 = 0
+        self.B16 = 0
+        self.B22 = 0
+        self.B26 = 0
+        self.B66 = 0
+
+
+    cpdef void make_smeared(Laminate self):
+        r"""Make a laminated with smeared properties
+
+        The `B_{ij}` terms of the constitutive matrix are set to zero.
+
+        The `D_{ij}` terms are calculated from the membrane terms `A_{ij}`
+        according to `D_{ij} = (h^2 A_{ij})/12`, where `h` is the
+        laminate thickness.
+
+        """
+        if self.offset != 0.:
+            raise NotImplementedError(
+                    'Laminates with offset cannot be made smeared!')
+
         self.B11 = 0
         self.B12 = 0
         self.B16 = 0
@@ -662,6 +687,14 @@ cdef class Laminate:
         self.B26 = 0
         self.B66 = 0
 
+        self.D11 = self.h**2/12 * self.A11
+        self.D12 = self.h**2/12 * self.A12
+        self.D16 = self.h**2/12 * self.A16
+        self.D22 = self.h**2/12 * self.A22
+        self.D26 = self.h**2/12 * self.A26
+        self.D66 = self.h**2/12 * self.A66
+
+
     cpdef LaminationParameters calc_lamination_parameters(Laminate self):
         r"""Calculate the lamination parameters.
 
@@ -715,11 +748,11 @@ cdef class Laminate:
         return lp
 
 
-cpdef LaminationParameters force_balanced_LP(LaminationParameters lp):
-    r"""Force balanced lamination parameters
+cpdef LaminationParameters make_balanced_LP(LaminationParameters lp):
+    r"""Make balanced lamination parameters
 
     The lamination parameters `\xi_{A2}` and `\xi_{A4}` are set to null to
-    force a balanced laminate.
+    make a balanced laminate.
 
     """
     lp.xiA2 = 0
@@ -727,10 +760,10 @@ cpdef LaminationParameters force_balanced_LP(LaminationParameters lp):
     return lp
 
 
-cpdef LaminationParameters force_symmetric_LP(LaminationParameters lp):
-    r"""Force symmetric lamination parameters
+cpdef LaminationParameters make_symmetric_LP(LaminationParameters lp):
+    r"""Make symmetric lamination parameters
 
-    The lamination parameters `\xi_{Bi}` are set to null to force a symmetric
+    The lamination parameters `\xi_{Bi}` are set to null to make a symmetric
     laminate.
 
     """
@@ -741,11 +774,11 @@ cpdef LaminationParameters force_symmetric_LP(LaminationParameters lp):
     return lp
 
 
-cpdef LaminationParameters force_orthotropic_LP(LaminationParameters lp):
-    r"""Force orthotropic lamination parameters
+cpdef LaminationParameters make_orthotropic_LP(LaminationParameters lp):
+    r"""Make orthotropic lamination parameters
 
     The lamination parameters `\xi_{A2}`, `\xi_{A4}`, `\xi_{B2}`, `\xi_{B4}`,
-    `\xi_{D2}` and `\xi_{D4}` are set to null to force an orthotropic laminate.
+    `\xi_{D2}` and `\xi_{D4}` are set to null to make an orthotropic laminate.
     The `\xi_{D2}` and `\xi_{D4}` are related to the bend-twist coupling and
     become often very small for balanced laminates with a large amount of
     plies.
diff --git a/composites/kassapoglou.py b/composites/kassapoglou.py
new file mode 100644
index 0000000..38ec4ce
--- /dev/null
+++ b/composites/kassapoglou.py
@@ -0,0 +1,218 @@
+r"""
+=====================================================
+Kassapoglou's methods (:mod:`composites.kassapoglou`)
+=====================================================
+
+.. currentmodule::composites.kassapoglou
+
+
+Methods based on Kassapoglou's book.
+
+
+Reference:
+
+    Kassapoglou. Design and Analysis of Composite Structures. 2nd Edition. John Wiley & Sons Ltd, 2013.
+
+
+"""
+import numpy as np
+from numpy import pi, tan, sqrt
+
+
+def calc_Nxx_crit(a, b, m, n, D11, D12, D22, D66):
+    r"""Calculate uniaxial compression buckling for a composite plate
+
+    The output of this function is the result of Eq. 6.6, section 6.2 page 129.
+    If `m` or `n` is set to ``None``, the function searchers for the critical number of
+    half-waves in the corresponding direction, up to 10 half-waves.
+
+    Reference:
+
+        Kassapoglou. Design and Analysis of Composite Structures. 2nd Edition. John Wiley & Sons Ltd, 2013.
+
+    Parameters
+    ----------
+    a, b : float
+        Plate length and width.
+    m, n : int or None
+        Number of half-waves along the plate length and width, respectively.
+    D11, D12, D22, D66 : float
+        Terms of the D matrix.
+
+    Result
+    ------
+    Nxx_crit : float
+        Critical uniaxial compression buckling load.
+
+    """
+    if m is None and n is None:
+        raise NotImplementedError("Only m or n can be None, not both")
+    if m is None:
+        m = np.arange(1, 11)
+    if n is None:
+        n = np.arange(1, 11)
+    AR = a/b
+    m = np.atleast_1d(m)
+    n = np.atleast_1d(n)
+    Nxx_crit = pi**2*(D11*m**4
+                      + 2*(D12 + 2*D66)*m**2*n**2*AR**2
+                      + D22*n**4*AR**4)/(a**2*m**2)
+    return Nxx_crit.min()
+
+
+def calc_Nxy_crit(a, D11, D12, D16, D22, D66, rtol=1e-5, atol=1e-6, max_iter=50):
+    r"""Calculate shear buckling for a composite plate
+
+    The output of this function is the result of Eq. 6.28, section 6.4 page
+    137. Variables `AR` and `\alpha` are solved using a Newton-Raphson scheme
+    that finds the solution of Eqs. 6.29 and 6.30 simultaneously.
+
+    Reference:
+
+        Kassapoglou. Design and Analysis of Composite Structures. 2nd Edition. John Wiley & Sons Ltd, 2013.
+
+    Parameters
+    ----------
+    a : float
+        Plate length.
+    D11, D12, D16, D22, D66 : float
+        Terms of the D matrix.
+    rtol, atol : float
+        Relative and absolute tolerances used to solve Eq. 6.30.
+    max_iter : int
+        Maximum number of iterations used in the Newton-Raphson scheme.
+
+    Result
+    ------
+    Nxy_crit : float
+        Critical shear buckling load.
+
+    """
+    alpha = np.pi/6
+    for i in range(max_iter):
+        AR = (D11/(D11*tan(alpha)**4 + 2*(D12 + 2*D66)*tan(alpha)**2 + D22))**(1/4)
+        expr = 3*D11*AR**4*tan(alpha)**4 + (6*D11*AR**2 + 2*(D12 + 2*D66)*AR**4)*tan(alpha)**2 - (D11 + 2*(D12 + 2*D66)*AR**2 + D22*AR**4)
+        if np.isclose(expr, 0, atol=1e-6, rtol=1e-5):
+            break
+        dexpr_dalpha = 3*AR**4*D11*(4*tan(alpha)**2 + 4)*tan(alpha)**3 + (AR**4*(2*D12 + 4*D66) + 6*AR**2*D11)*(2*tan(alpha)**2 + 2)*tan(alpha)
+        dalpha = -expr/dexpr_dalpha
+        alpha += dalpha
+    else:
+        raise RuntimeError('Newton-Raphson scheme did not converge!')
+    Nxy_crit = pi**2/(2*AR**2*a**2*tan(alpha))*(
+        D11*(1 + 6*tan(alpha)**2*AR**2 + tan(alpha)**4*AR**4)
+        + 2*(D12 + 2*D66)*(AR**2 + AR**4*tan(alpha)**2)
+        + D22*AR**4)
+    return Nxy_crit
+
+
+def calc_beff(b, Px, Pcr, A11, A12, A22):
+    r"""Calculate the effective width of a plate
+
+    The output of this function is the result of Eq. 7.15, section 7.1 page
+    164. The effective width `b_{eff}` is defined as:
+
+    `\int_y N_{xx} dy = 2{N_x}_{max} b_{eff}`
+
+
+    Reference:
+
+        Kassapoglou. Design and Analysis of Composite Structures. 2nd Edition. John Wiley & Sons Ltd, 2013.
+
+    Parameters
+    ----------
+    b : float
+        Plate width.
+    Px : float
+        Applied compressive force.
+    Pcr : float
+        Critical buckling force.
+    A11, A12, A22: float
+        Terms of the A matrix.
+
+    Result
+    ------
+    beff : float
+        Effective width of the plate.
+
+    """
+    beff = b/(2*(1 + 2*(1 + A12/A11)*(1 - Pcr/Px)*A11/(A11 + 3*A22)))
+    return beff
+
+
+def calc_Nxx_crit_combined_shear(k, a, b, D11, D12, D22, D66):
+    r"""Calculate combined uniaxial-shear buckling load for a composite plate
+
+    The output of this function is the result of Eq. 6.34, section 6.5 page
+    142. This function calculates the critical `N_{xx}` buckling load under the
+    current level of shear load given by `N_{xy} = k N_{xx}`.
+
+    Reference:
+
+        Kassapoglou. Design and Analysis of Composite Structures. 2nd Edition. John Wiley & Sons Ltd, 2013.
+
+    Parameters
+    ----------
+    k : float
+        Load ratio defined as `k=N_{xy}/N_{xx}`, where `N_{xy}` is the current
+        shear load, and `N_{xx}` the current compression load.
+    a, b : float
+        Plate length and width.
+    D11, D12, D22, D66 : float
+        Terms of the D matrix.
+
+    Result
+    ------
+    N0 : float
+        Critical `N_{xx}` buckling load under the current level of shear load
+        given by `N_{xy} = k N_{xx}`.
+
+    """
+    den = 2 - 8192/81*a**2*k**2/(b**2*pi**4)
+    term = 9 + 65536/81*a**2*k**2/(b**2*pi**4)
+    rhs_term1 = 5 + sqrt(term)
+    rhs_term2 = 5 - sqrt(term)
+    Nxx_crit1 = pi**2/a**2*(D11 + 2*(D12 + 2*D66)*a**2/b**2 + D22*a**4/b**4)/den*rhs_term1
+    Nxx_crit2 = pi**2/a**2*(D11 + 2*(D12 + 2*D66)*a**2/b**2 + D22*a**4/b**4)/den*rhs_term2
+    return min(abs(Nxx_crit1), abs(Nxx_crit2))
+
+
+def calc_Nxx_crit_combined_shear_full(Nxy, a, b, D11, D12, D16, D22, D26, D66):
+    r"""Calculate combined uniaxial-shear buckling load for a composite plate
+
+    This solution does not ignore the D16 and D26 terms. Furthermore, this
+    solution takes as input `N_{xy}` instead of the `k=N_{xy}/N_{xx}`
+    originally used by Kassapoglou.
+
+    Based on section 6.5. This function calculates the critical `N_{xx}` by
+    `N_{xy}`.
+
+    Reference:
+
+        Kassapoglou. Design and Analysis of Composite Structures. 2nd Edition. John Wiley & Sons Ltd, 2013.
+
+    Parameters
+    ----------
+    Nxy : float
+        Shear load `N_{xy}`.
+    a, b : float
+        Plate length and width.
+    D11, D12, D16, D22, D26, D66 : float
+        All terms of the D matrix.
+
+    Result
+    ------
+    Nxx_crit : float
+        Critical `N_{xx}` buckling load under the current level of shear load
+        given by `N_{xy}`.
+
+    """
+    a11 = pi**4*D11*b/(4*a**3) + pi**4*D12/(2*a*b) + pi**4*D22*a/(4*b**3) + pi**4*D66/(a*b)
+    a12 = -160*pi**2*D16/(9*a**2) - 160*pi**2*D26/(9*b**2)
+    a21 = -160*pi**2*D16/(9*a**2) - 160*pi**2*D26/(9*b**2)
+    a22 = 4*pi**4*D11*b/a**3 + 8*pi**4*D12/(a*b) + 4*pi**4*D22*a/b**3 + 16*pi**4*D66/(a*b)
+    term = 16384*Nxy**2 + 4608*Nxy*a12 + 4608*Nxy*a21 + 1296*a11**2 - 648*a11*a22 + 1296*a12*a21 + 81*a22**2
+    Nxx_crit1 = a*(36*a11 + 9*a22 - sqrt(term))/(18*pi**2*b)
+    Nxx_crit2 = a*(36*a11 + 9*a22 + sqrt(term))/(18*pi**2*b)
+    return min(abs(Nxx_crit1), abs(Nxx_crit2))
+
diff --git a/doc/source/api.rst b/doc/source/api.rst
new file mode 100644
index 0000000..d6fbe67
--- /dev/null
+++ b/doc/source/api.rst
@@ -0,0 +1,2 @@
+.. automodule:: composites
+    :members:
diff --git a/doc/source/index.rst b/doc/source/index.rst
index edf499b..d95ff5d 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -1,8 +1,8 @@
 Documentation for the ``composites`` module
 ===========================================
 
-High-performance module to calculate properties of laminated composite
-materials. Usually, this module is used to calculate:
+High-performance module to calculate properties and perform analysis on
+composites. With the ``composites`` module, you are able to calculate:
 
 * A, B, D, E plane-stress stiffness matrices for plates
     - A, B, D, for classical plate theory (CLT, or CLPT)
@@ -14,6 +14,9 @@ materials. Usually, this module is used to calculate:
 
 * Stiffness matrices (ABDE) based on lamination parameters
 
+* Based on Kassapoglou's book, local buckling under compression, shear, and
+  post-buckling 
+
 
 Code repository
 ---------------
@@ -25,22 +28,24 @@ Citing this library
 -------------------
 
 Castro, S. G. P. Methods to calculate composite plate properties (Version
-0.6.5) [Computer software]. 2024. https://doi.org/10.5281/zenodo.2871782
+0.7.0) [Computer software]. 2024. https://doi.org/10.5281/zenodo.2871782
 
 Bibtex :
     
     @misc{composites2024,
         author = {Castro, Saullo G. P.},
         doi = {10.5281/zenodo.2871782},
-        title = {{Methods to calculate composite plate properties (Version 0.6.5) [Computer software]. 2024}}
+        title = {{Methods to calculate composite plate properties (Version 0.7.0) [Computer software]. 2024}}
         }
 
 
-Classes and functions
----------------------
+composites API
+--------------
+
+.. toctree::
+    :maxdepth: 1
 
-.. automodule:: composites
-    :members:
+    api.rst
 
 
 License
diff --git a/tests/test_core.py b/tests/test_core.py
index 742e8bb..b053ed1 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -6,9 +6,10 @@
 from composites.utils import (read_laminaprop, laminated_plate,
         isotropic_plate)
 from composites.core import (laminate_from_LaminationParameters,
-        laminate_from_lamination_parameters, force_balanced_LP,
-        force_orthotropic_LP, force_symmetric_LP, Lamina,
-        GradABDE, LaminationParameters)
+                             laminate_from_lamination_parameters,
+                             make_balanced_LP, make_orthotropic_LP,
+                             make_symmetric_LP, Lamina, GradABDE,
+                             LaminationParameters)
 
 
 def test_lampar_tri_axial():
@@ -47,7 +48,7 @@ def test_lampar_tri_axial():
     E = np.array([[2.66917293e+10,  0.00000000e+00],
                   [0.00000000e+00,  2.66917293e+10]])
     assert np.allclose(lam.A, A)
-    lam.force_symmetric()
+    lam.make_symmetric()
     assert np.allclose(lam.B, B)
     assert np.allclose(lam.D, D)
     assert np.allclose(lam.E, E)
@@ -94,7 +95,7 @@ def test_lampar_plane_stress():
     E = np.array([[2.66917293e+10, 0.00000000e+00],
                   [0.00000000e+00, 2.66917293e+10]])
     assert np.allclose(lam.A, A)
-    lam.force_symmetric()
+    lam.make_symmetric()
     assert np.allclose(lam.B, B)
     assert np.allclose(lam.D, D)
     assert np.allclose(lam.E, E)
@@ -159,16 +160,16 @@ def test_laminated_plate_plane_stress():
     assert np.allclose(lam_2.D, lam.D)
     assert np.allclose(lam_2.E, lam.E)
 
-    lam.force_balanced()
-    force_balanced_LP(lp)
+    lam.make_balanced()
+    make_balanced_LP(lp)
     lam_2 = laminate_from_LaminationParameters(thickness, matlamina, lp)
     assert np.allclose(lam_2.A, lam.A)
     assert np.allclose(lam_2.B, lam.B)
     assert np.allclose(lam_2.D, lam.D)
     assert np.allclose(lam_2.E, lam.E)
 
-    lam.force_orthotropic()
-    force_orthotropic_LP(lp)
+    lam.make_orthotropic()
+    make_orthotropic_LP(lp)
     lam_2 = laminate_from_LaminationParameters(thickness, matlamina, lp)
     assert np.allclose(lam_2.A, lam.A)
     assert np.allclose(lam_2.B, lam.B)
@@ -177,14 +178,20 @@ def test_laminated_plate_plane_stress():
 
     lam = laminated_plate(stack, plyt, lamprop)
     lp = lam.calc_lamination_parameters()
-    lam.force_symmetric()
-    force_symmetric_LP(lp)
+    lam.make_symmetric()
+    make_symmetric_LP(lp)
     lam_2 = laminate_from_LaminationParameters(thickness, matlamina, lp)
     assert np.allclose(lam_2.A, lam.A)
     assert np.allclose(lam_2.B, lam.B)
     assert np.allclose(lam_2.D, lam.D)
     assert np.allclose(lam_2.E, lam.E)
 
+    lam = laminated_plate(stack, plyt, lamprop)
+    lp = lam.calc_lamination_parameters()
+    lam.make_smeared()
+    assert np.allclose(lam.D, lam.h**2/12*lam.A)
+    assert np.allclose(lam.B, 0)
+
 
 def test_isotropic_plate():
     E = 71e9
@@ -214,15 +221,19 @@ def test_errors():
     lam = isotropic_plate(thickness=thick, E=E, nu=nu)
     lam.offset = 1.
     try:
-        lam.force_balanced()
+        lam.make_balanced()
+    except RuntimeError:
+        pass
+    try:
+        lam.make_orthotropic()
     except RuntimeError:
         pass
     try:
-        lam.force_orthotropic()
+        lam.make_symmetric()
     except RuntimeError:
         pass
     try:
-        lam.force_symmetric()
+        lam.make_smeared()
     except RuntimeError:
         pass
     try:
diff --git a/tests/test_kassapoglou.py b/tests/test_kassapoglou.py
new file mode 100644
index 0000000..c700106
--- /dev/null
+++ b/tests/test_kassapoglou.py
@@ -0,0 +1,192 @@
+import sys
+sys.path.append('..')
+
+import numpy as np
+from numpy import pi, sqrt
+
+from composites import laminated_plate
+from composites.kassapoglou import (calc_Nxx_crit,
+                                    calc_Nxy_crit,
+                                    calc_Nxx_crit_combined_shear,
+                                    calc_Nxx_crit_combined_shear_full,
+                                    calc_beff,
+                                    )
+
+
+def test_calc_Nxx_crit():
+    """Verificatoin based on Kassapoglou's Figs. 6.4 and 6.5
+
+    """
+    alu = laminated_plate(stack=[0], plyt=0.5715e-3, laminaprop=(70e9, 0.3))
+    D11, D12, D22, D66 = 0.66, 0.47, 0.66, 0.49
+    AR_values = np.linspace(0.6, 2.8)
+    res1 = []
+    res2 = []
+    res3 = []
+    res_alu1 = []
+    res_alu2 = []
+    res_alu3 = []
+    b = 0.508/2
+    for AR in AR_values:
+        a = AR*b
+        res1.append(calc_Nxx_crit(a, b, 1, 1, D11, D12, D22, D66))
+        res2.append(calc_Nxx_crit(a, b, 2, 1, D11, D12, D22, D66))
+        res3.append(calc_Nxx_crit(a, b, 3, 1, D11, D12, D22, D66))
+
+        res_alu1.append(calc_Nxx_crit(a, b, 1, 1, alu.D11, alu.D12, alu.D22, alu.D66))
+        res_alu2.append(calc_Nxx_crit(a, b, 2, 1, alu.D11, alu.D12, alu.D22, alu.D66))
+        res_alu3.append(calc_Nxx_crit(a, b, 3, 1, alu.D11, alu.D12, alu.D22, alu.D66))
+
+    res1_expected = [760.4492762134736, 728.400368245419, 703.8761399336366, 685.191127113606, 671.1307610092099, 660.8018833650219, 653.5358314008664, 648.8239785928772, 646.2738786053692, 645.5788159767615, 646.4962762781008, 648.8324700409954, 652.4310404371092, 657.1647103179126, 662.9290256450245, 669.6376148748768, 677.2185585936345, 685.6115818818151, 694.7658630211979, 704.6383086268357, 715.192185097175, 726.3960246767871, 738.2227449149805, 750.6489352387343, 763.6542753519922, 777.2210583401934, 791.3337974784167, 805.9789003643929, 821.1443975169391, 836.8197152790613, 852.995484948758, 869.6633816801207, 886.8159879639275, 904.4466774933751, 922.5495160089763, 941.119176343723, 960.1508653910068, 979.6402611206446, 999.5834580935449, 1019.9769201892686, 1040.817439475499, 1062.102100324078, 1083.828248022461, 1105.9934612482962, 1128.5955278731656, 1151.6324236431244, 1175.1022933516867, 1199.003434177715, 1223.3342809083138, 1248.093392806889]
+    res2_expected = [1574.5738453392305, 1425.2156696076072, 1304.4297339276068, 1205.4741816411472, 1123.4907376437611, 1054.9067689137432, 997.0476243303884, 947.8787977982346, 905.8305039745393, 869.67588101298, 838.4448711977442, 811.3623166552635, 787.8027900721946, 767.2571828544187, 749.3076788484112, 733.6087918799004, 719.8728442935457, 707.8587364114169, 697.3631813606328, 688.2138056014026, 680.2636747275131, 673.3869177172485, 667.4752047678446, 662.4348935872179, 658.1847029911407, 654.6538053213732, 651.7802536782276, 649.5096784526276, 647.7942017198433, 646.5915288518977, 645.8641850407839, 645.5788709028703, 645.7059164012659, 646.2188163087604, 647.0938335874052, 648.3096595691644, 649.8471218276071, 651.6889322420013, 653.819469055981, 656.224587787787, 658.8914567081547, 661.8084133044526, 664.9648387265003, 668.3510476848932, 671.9581916659561, 675.7781736539106, 679.8035728228145, 684.0275778881181, 688.4439279982369, 693.0468602067966]
+    res3_expected = [2971.8346291977537, 2633.2313905592655, 2358.731949245679, 2133.1671277763076, 1945.6058071661016, 1788.0095631162233, 1654.360430605087, 1540.0807709179026, 1441.638567033499, 1356.2733797976512, 1281.8025793497804, 1216.4820604743586, 1158.9046117163628, 1107.9247387390258, 1062.602355196917, 1022.1601171979253, 985.9507500175343, 953.4317793804016, 924.1458098225442, 897.7050008794315, 873.7787501370861, 852.0838477966455, 832.3765518018913, 814.446166995635, 798.1098107133652, 783.2081207240418, 769.6017165033668, 757.1682664309185, 745.8000451760562, 735.4018898258297, 725.8894820622249, 717.1878982720699, 709.2303808723597, 701.9572931018719, 695.315226625222, 689.2562369393289, 683.737186084723, 678.7191757897579, 674.1670571026093, 670.0490049393197, 666.3361479090431, 663.0022453583597, 660.0234048743622, 657.3778345558914, 655.0456252471773, 653.0085586627191, 651.2499379441467, 649.7544377012231, 648.5079710178828, 647.4975712647844]
+
+    assert np.allclose(res1, res1_expected)
+    assert np.allclose(res2, res2_expected)
+    assert np.allclose(res3, res3_expected)
+
+    #import matplotlib.pyplot as plt
+    #plt.plot(AR_values, np.asarray(res1)/1000, label='m=1')
+    #plt.plot(AR_values, np.asarray(res2)/1000, label='m=2')
+    #plt.plot(AR_values, np.asarray(res3)/1000, label='m=3')
+    #plt.plot(AR_values, np.asarray(res_alu1)/1000, label='m=1')
+    #plt.plot(AR_values, np.asarray(res_alu2)/1000, label='m=2')
+    #plt.plot(AR_values, np.asarray(res_alu3)/1000, label='m=3')
+    #plt.ylim(-0.1, 0.9)
+    #plt.legend()
+    #plt.xlabel('Length/Width')
+    #plt.ylabel('Buckling load ${N_{xx}}_{crit}$ [N/mm]')
+    #plt.title("Reproducing Kassapoglou's Figs. 6.4 and 6.5")
+    #plt.show()
+
+    AR_values = np.linspace(0.6, 5, 200)
+    res = [calc_Nxx_crit(b*AR, b, None, 1, D11, D12, D22, D66) for AR in
+           AR_values]
+    #import matplotlib.pyplot as plt
+    #plt.plot(AR_values, np.asarray(res))
+    #plt.show()
+
+    assert np.allclose(calc_Nxx_crit(b*AR, b, 1, 1, D11, D12, D22, D66),
+                       calc_Nxx_crit(b*AR, b, 1, None, D11, D12, D22, D66))
+
+    try:
+        calc_Nxx_crit(b*AR, b, None, None, D11, D12, D22, D66)
+    except NotImplementedError:
+        pass
+
+
+def test_calc_Nxy_crit():
+    """Verificatoin based on Kassapoglou's Fig. 6.11
+
+    """
+    laminaprop = (68.9e9, 68.0e9, 0.05, 4.83e9, 4.83e9, 4.83e9)
+    lam = laminated_plate(stack=[0, 90]*8, plyt=0.1905e-3/2, laminaprop=laminaprop)
+    a_values = np.linspace(0.1, 0.5)
+    res = []
+    for a in a_values:
+        res.append(calc_Nxy_crit(a, lam.D11, lam.D12, lam.D16, lam.D22, lam.D66))
+    res_expected = [79688.55992406543, 68114.0022704454, 58889.57598574364, 51419.57333450177, 45285.73547400736, 40187.4044061502, 35903.96554281875, 32270.57385354717, 29162.053403091162, 26481.969879263816, 24155.06026735022, 22121.89066686104, 20335.02310316517, 18756.22315240477, 17354.39749457425, 16104.051205932257, 14984.120320908534, 13977.078850002274, 13068.248915899261, 12245.262872171586, 11497.640308736318, 10816.452732075362, 10194.055750316005, 9623.87366720392, 9100.225083361767, 8618.180819678446, 8173.447493599944, 7762.271588205649, 7381.359992966362, 7027.813861439158, 6699.073294971503, 6392.870873657024, 6107.192453563187, 5840.243960125795, 5590.423152014058, 5356.295523016742, 5136.573663123334, 4930.099522731353, 4735.829122489075, 4552.8193308193, 4380.216395633825, 4217.245969223063, 4063.204408198967, 3917.45116557157, 3779.4021210406136, 3648.5237195644836, 3524.327808168894, 3406.3670775281935, 3294.2310286958063, 3187.5423969626177]
+    assert np.allclose(res, res_expected)
+
+    try:
+        calc_Nxy_crit(a, -lam.D11, lam.D12, lam.D16, lam.D22, lam.D66)
+    except RuntimeError:
+        pass
+
+    #import matplotlib.pyplot as plt
+    #plt.plot(a_values*1000, np.asarray(res)/1000)
+    #plt.xlabel('Plate dimension a [mm]')
+    #plt.ylabel('Buckling load [N/mm]')
+    #plt.title("Reproducing Kassapoglou's Fig. 6.11")
+
+
+def test_calc_Nxx_crit_combined_shear_full():
+    D11, D12, D22, D66 = 0.66, 0.47, 0.66, 0.49
+    b = 0.508/2
+    a = 2*b
+    Nxx_crit = calc_Nxx_crit(a, b, 1, 1, D11, D12, D22, D66)
+    Nxy_crit = calc_Nxy_crit(a, D11, D12, 0, D22, D66)
+    x_values = []
+    y_values = []
+    for Nxy in np.linspace(0.01*Nxy_crit, 0.99*Nxy_crit):
+        Nxx = calc_Nxx_crit_combined_shear_full(Nxy, a, b, D11, D12, 0, D22, 0, D66)
+        assert np.isclose(Nxx, calc_Nxx_crit_combined_shear(Nxy/Nxx, a, b, D11, D12, D22, D66))
+        x_values.append(Nxx/Nxx_crit)
+        y_values.append(Nxy/Nxy_crit)
+
+    #import matplotlib.pyplot as plt
+    #plt.plot(x_values, y_values)
+    #plt.ylim(0., 1.1)
+    #plt.xlabel('Length/Width')
+    #plt.ylabel('Buckling load ${N_{xx}}_{crit}$ [N/mm]')
+    #plt.title("Reproducing Kassapoglou's Figs. 6.4 and 6.5")
+    #plt.show()
+
+
+def test_calc_beff():
+    """Verificatoin based on Kassapoglou's Fig. 7.10
+
+    """
+    laminaprop = (68.9e9, 68.0e9, 0.05, 4.83e9, 4.83e9, 4.83e9)
+    plyt = 0.1905e-3/2/16
+    lam = laminated_plate(stack=[0]*8, plyt=plyt, laminaprop=laminaprop)
+
+    Px = np.linspace(10, 1, 40)
+    Pcr = 1
+
+    a = b = 1
+
+    A11 = lam.A11
+    A22 = A11*10.
+    A12 = A11*0.7
+    beffs1 = calc_beff(b, Px, Pcr, A11, A12, A22)
+
+    A11 = lam.A11
+    A22 = A11*1.
+    A12 = A11*0.3
+    beffs2 = calc_beff(b, Px, Pcr, A11, A12, A22)
+
+    A11 = lam.A11
+    A22 = A11*0.1
+    A12 = A11*0.05
+    beffs3 = calc_beff(b, Px, Pcr, A11, A12, A22)
+
+    beffs1_expected = [0.45507927187316505, 0.4551866068538131, 0.45529918984223244, 0.4554174153756495, 0.45554171855541725, 0.4556725803977709, 0.4558105340545082, 0.4559561720731905, 0.4561101549053356, 0.45627322092026296, 0.45644619824489074, 0.45663001883010385, 0.45682573524800907, 0.45703454085930917, 0.457257794166946, 0.45749704840613936, 0.45775408773319887, 0.45803097179621455, 0.45833009104349853, 0.4586542359140656, 0.45900668415263557, 0.45939131204105105, 0.4598127375544637, 0.460276506669263, 0.4607893388006151, 0.4613594544864163, 0.46199701937406856, 0.4627147556963568, 0.4635287998959824, 0.46445993031358884, 0.4655353656705211, 0.4667914699658148, 0.4682779456193354, 0.4700645666210135, 0.47225244831338414, 0.4749938710468252, 0.47852932921695196, 0.48326222513948147, 0.48992493085736866, 0.5]
+    beffs2_expected = [0.31545741324921134, 0.31576330183988066, 0.31608462911037466, 0.31642259414225943, 0.3167785234899329, 0.31715388858246, 0.31755032605613837, 0.31796966161026835, 0.31841393811955543, 0.3188854489164087, 0.31938677738741617, 0.31992084432717677, 0.320490964882373, 0.3211009174311927, 0.3217550274223035, 0.3224582701062215, 0.3232163973196689, 0.3240360951599672, 0.32492518170158186, 0.32589285714285715, 0.3269500233535731, 0.3281096963761019, 0.3293875450334534, 0.33080260303687636, 0.332378223495702, 0.33414337788578374, 0.3361344537815126, 0.3383977900552486, 0.34099332839140106, 0.34400000000000003, 0.3475238922675934, 0.35171102661596954, 0.3567681007345226, 0.3629976580796253, 0.37086092715231794, 0.38109756097560976, 0.3949730700179534, 0.4148471615720524, 0.445682451253482, 0.5]
+    beffs3_expected = [0.2037617554858934, 0.20407911001236095, 0.20441288359117424, 0.20476438427492838, 0.20513506285102973, 0.20552653285675004, 0.20594059405940593, 0.20637926012234198, 0.20684479135394776, 0.20733973366367295, 0.20786696514230893, 0.20842975206611566, 0.20903181662675333, 0.20967741935483866, 0.2103714600956234, 0.2111196025983951, 0.21192842942345927, 0.21280563613758807, 0.21376027693639116, 0.2148030783159801, 0.21594684385382062, 0.21720698254364088, 0.2186022070415134, 0.22015546918378676, 0.22189523248969983, 0.22385723231058233, 0.22608695652173913, 0.22864321608040203, 0.23160340821068942, 0.23507148864592092, 0.23919043238270468, 0.24416243654822334, 0.25028312570781425, 0.25800256081946227, 0.2680412371134021, 0.2816291161178509, 0.30105263157894746, 0.3310991957104557, 0.3837638376383766, 0.5]
+
+    assert np.allclose(beffs1, beffs1_expected)
+    assert np.allclose(beffs2, beffs2_expected)
+    assert np.allclose(beffs3, beffs3_expected)
+
+    #import matplotlib.pyplot as plt
+    #plt.plot(Pcr/Px, beffs1/b)
+    #plt.plot(Pcr/Px, beffs2/b)
+    #plt.plot(Pcr/Px, beffs3/b)
+    #plt.show()
+
+
+if __name__ == '__main__':
+    test_calc_Nxx_crit()
+    test_calc_Nxy_crit()
+    test_calc_Nxx_crit_combined_shear_full()
+    test_calc_beff()
+
+    if False:
+        laminaprop = (68.9e9, 68.0e9, 0.05, 4.83e9, 4.83e9, 4.83e9)
+        plyt = 0.1905e-3/2/16
+        lam = laminated_plate(stack=[0]*8, plyt=plyt, laminaprop=laminaprop)
+
+        Px = np.linspace(10, 1, 40)
+        Pcr = 1
+
+        b = 1
+
+        A11 = lam.A11
+        A22 = A11*10.
+        A12 = A11*0.7
+        beffs1 = calc_beff(b, Px, Pcr, A11, A12, A22)
+
+        import matplotlib.pyplot as plt
+        plt.plot(Pcr/Px, beffs1/a)
+        plt.show()
diff --git a/theory/kassapoglou.py b/theory/kassapoglou.py
new file mode 100644
index 0000000..34d9b5b
--- /dev/null
+++ b/theory/kassapoglou.py
@@ -0,0 +1,78 @@
+print("""
+# shear buckling
+# section 6.4
+# page 138, Eq. 6.30
+# differentiation used in the Newton-Raphson approach herein implemented
+""")
+import sympy
+from sympy import pi, sin, tan
+
+sympy.var("D11, D12, D22, D66, AR, alpha")
+expr = (3*D11*AR**4*tan(alpha)**4
+        + (6*D11*AR**2 + 2*(D12 + 2*D66)*AR**4)*tan(alpha)**2
+        - (D11 + 2*(D12 + 2*D66)*AR**2 + D22*AR**4))
+dexpr_dalpha = expr.diff(alpha)
+print("dexpr_dalpha =", dexpr_dalpha)
+print()
+print("""
+# combined compression-shear buckling
+# section 6.5
+# Ritz method therein used
+""")
+sympy.var("D11, D12, D16, D22, D26, D66")
+sympy.var("Nxx, Nxy, k")
+#Nxy = k*Nxx
+sympy.var("w1, w2")
+sympy.var("a, b, x, y", positive=True)
+w = w1*sin(pi*x/a)*sin(pi*y/b) + w2*sin(2*pi*x/a)*sin(2*pi*y/b)
+#Eq. 6.32
+integrand_Uc = (D11*(w.diff(x, x))**2
+            + 2*D12*w.diff(x, x)*w.diff(y, y)
+              + D22*w.diff(y, y)**2
+            + 4*D66*w.diff(x, y)**2
+            + 4*D16*w.diff(x, x)*w.diff(x, y)
+            + 4*D26*w.diff(y, y)*w.diff(x, y)
+               )/2
+integrand_Vc = (
+                Nxx*w.diff(x)**2/2
+              + Nxy*w.diff(x)*w.diff(y)
+               )
+PIc = sympy.integrate(integrand_Uc + integrand_Vc, (x, 0, a), (y, 0, b))
+eq1 = PIc.diff(w1)
+eq2 = PIc.diff(w2)
+#sympy.solve((eq1, eq2), (w1, w2))
+#
+# eq1.collect((w1, w2))
+# eq1, terms for w1
+print(eq1.subs(dict(w2=0)).collect((w1,Nxx,Nxy)))
+print(eq1.subs(dict(w1=0)).collect((w2,Nxx,Nxy)))
+a11 = pi**4*D11*b/(4*a**3) + pi**4*D12/(2*a*b) + pi**4*D22*a/(4*b**3) + pi**4*D66/(a*b)
+a12 = -160*pi**2*D16/(9*a**2) - 160*pi**2*D26/(9*b**2)
+b11 = pi**2*Nxx*b/(4*a)
+b12 = -32*Nxy/9
+
+print(eq2.subs(dict(w2=0)).collect((w1,Nxx,Nxy)))
+print(eq2.subs(dict(w1=0)).collect((w2,Nxx,Nxy)))
+a21 = -160*pi**2*D16/(9*a**2) - 160*pi**2*D26/(9*b**2)
+a22 = 4*pi**4*D11*b/a**3 + 8*pi**4*D12/(a*b) + 4*pi**4*D22*a/b**3 + 16*pi**4*D66/(a*b)
+b21 = -32*Nxy/9
+b22 = pi**2*Nxx*b/a
+
+print('a11 =', a11)
+print('a12 =', a12)
+print('a21 =', a21)
+print('a22 =', a22)
+#
+a11, a12, a21, a22 = sympy.var('a11, a12, a21, a22')
+A = sympy.Matrix([[a11, a12],
+                  [a21, a22]])
+B = sympy.Matrix([[b11, b12],
+                  [b21, b22]])
+cp = sympy.det(A - B)
+N0 = sympy.solve(cp, Nxx)
+print('________')
+print()
+print('Solution')
+print('________')
+print(N0[0].expand().simplify())
+print(N0[1].expand().simplify())