def calc_scf(self):
        """Calculate improved shear correction factors

        Reference:

            Vlachoutsis, S. "Shear correction factors for plates and shells",
            Int. Journal for Numerical Methods in Engineering, Vol. 33,
            1537-1552, 1992.

            http://onlinelibrary.wiley.com/doi/10.1002/nme.1620330712/full


        Using "one shear correction factor" (see reference), assuming:

        - constant G13, G23, E1, E2, nu12, nu21 within each ply
        - g1 calculated using z at the middle of each ply
        - zn1 = Laminate.offset

        Returns
        -------

        k13, k23 : tuple
            Shear correction factors. Also updates attributes: `scf_k13`
            and `scf_k23`.

        """
        D1 = 0
        R1 = 0
        den1 = 0

        D2 = 0
        R2 = 0
        den2 = 0

        offset = self.offset
        zbot = -self.h/2 + offset
        z1 = zbot

        for ply in self.plies:
            z2 = z1 + ply.h
            e1 = (ply.matobj.e1 * np.cos(np.deg2rad(ply.theta)) +
                  ply.matobj.e2 * np.sin(np.deg2rad(ply.theta)))
            e2 = (ply.matobj.e2 * np.cos(np.deg2rad(ply.theta)) +
                  ply.matobj.e1 * np.sin(np.deg2rad(ply.theta)))
            nu12 = (ply.matobj.nu12 * np.cos(np.deg2rad(ply.theta)) +
                  ply.matobj.nu21 * np.sin(np.deg2rad(ply.theta)))
            nu21 = (ply.matobj.nu21 * np.cos(np.deg2rad(ply.theta)) +
              ply.matobj.nu12 * np.sin(np.deg2rad(ply.theta)))

            D1 += e1 / (1 - nu12*nu21)
            R1 += D1*((z2 - offset)**3/3 - (z1 - offset)**3/3)
            g13 = ply.matobj.g13
            d1 = g13 * ply.h
            den1 += d1 * (self.h / ply.h) * D1**2*(15*offset*z1**4 + 30*offset*z1**2*zbot*(2*offset - zbot) - 15*offset*z2**4 + 30*offset*z2**2*zbot*(-2*offset + zbot) - 3*z1**5 + 10*z1**3*(-2*offset**2 - 2*offset*zbot + zbot**2) - 15*z1*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2) + 3*z2**5 + 10*z2**3*(2*offset**2 + 2*offset*zbot - zbot**2) + 15*z2*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2))/(60*g13)

            D2 += e2 / (1 - nu12*nu21)
            R2 += D2*((z2 - self.offset)**3/3 - (z1 - self.offset)**3/3)
            g23 = ply.matobj.g23
            d2 = g23 * ply.h
            den2 += d2 * (self.h / ply.h) * D2**2*(15*offset*z1**4 + 30*offset*z1**2*zbot*(2*offset - zbot) - 15*offset*z2**4 + 30*offset*z2**2*zbot*(-2*offset + zbot) - 3*z1**5 + 10*z1**3*(-2*offset**2 - 2*offset*zbot + zbot**2) - 15*z1*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2) + 3*z2**5 + 10*z2**3*(2*offset**2 + 2*offset*zbot - zbot**2) + 15*z2*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2))/(60*g23)

            z1 = z2

        self.scf_k13 = R1**2 / den1
        self.scf_k23 = R2**2 / den2

        return self.scf_k13, self.scf_k23
def calc_scf ( self ) : D1 = 0 R1 = 0 den1 = 0 D2 = 0 R2 = 0 den2 = 0 offset = self . offset zbot = - self . h / 2 + offset z1 = zbot for ply in self . plies : z2 = z1 + ply . h e1 = ( ply . matobj . e1 * np . cos ( np . deg2rad ( ply . theta ) ) + ply . matobj . e2 * np . sin ( np . deg2rad ( ply . theta ) ) ) e2 = ( ply . matobj . e2 * np . cos ( np . deg2rad ( ply . theta ) ) + ply . matobj . e1 * np . sin ( np . deg2rad ( ply . theta ) ) ) nu12 = ( ply . matobj . nu12 * np . cos ( np . deg2rad ( ply . theta ) ) + ply . matobj . nu21 * np . sin ( np . deg2rad ( ply . theta ) ) ) nu21 = ( ply . matobj . nu21 * np . cos ( np . deg2rad ( ply . theta ) ) + ply . matobj . nu12 * np . sin ( np . deg2rad ( ply . theta ) ) ) D1 += e1 / ( 1 - nu12 * nu21 ) R1 += D1 * ( ( z2 - offset ) ** 3 / 3 - ( z1 - offset ) ** 3 / 3 ) g13 = ply . matobj . g13 d1 = g13 * ply . h den1 += d1 * ( self . h / ply . h ) * D1 ** 2 * ( 15 * offset * z1 ** 4 + 30 * offset * z1 ** 2 * zbot * ( 2 * offset - zbot ) - 15 * offset * z2 ** 4 + 30 * offset * z2 ** 2 * zbot * ( - 2 * offset + zbot ) - 3 * z1 ** 5 + 10 * z1 ** 3 * ( - 2 * offset ** 2 - 2 * offset * zbot + zbot ** 2 ) - 15 * z1 * zbot ** 2 * ( 4 * offset ** 2 - 4 * offset * zbot + zbot ** 2 ) + 3 * z2 ** 5 + 10 * z2 ** 3 * ( 2 * offset ** 2 + 2 * offset * zbot - zbot ** 2 ) + 15 * z2 * zbot ** 2 * ( 4 * offset ** 2 - 4 * offset * zbot + zbot ** 2 ) ) / ( 60 * g13 ) D2 += e2 / ( 1 - nu12 * nu21 ) R2 += D2 * ( ( z2 - self . offset ) ** 3 / 3 - ( z1 - self . offset ) ** 3 / 3 ) g23 = ply . matobj . g23 d2 = g23 * ply . h den2 += d2 * ( self . h / ply . h ) * D2 ** 2 * ( 15 * offset * z1 ** 4 + 30 * offset * z1 ** 2 * zbot * ( 2 * offset - zbot ) - 15 * offset * z2 ** 4 + 30 * offset * z2 ** 2 * zbot * ( - 2 * offset + zbot ) - 3 * z1 ** 5 + 10 * z1 ** 3 * ( - 2 * offset ** 2 - 2 * offset * zbot + zbot ** 2 ) - 15 * z1 * zbot ** 2 * ( 4 * offset ** 2 - 4 * offset * zbot + zbot ** 2 ) + 3 * z2 ** 5 + 10 * z2 ** 3 * ( 2 * offset ** 2 + 2 * offset * zbot - zbot ** 2 ) + 15 * z2 * zbot ** 2 * ( 4 * offset ** 2 - 4 * offset * zbot + zbot ** 2 ) ) / ( 60 * g23 ) z1 = z2 self . scf_k13 = R1 ** 2 / den1 self . scf_k23 = R2 ** 2 / den2 return self . scf_k13 , self . scf_k23
Calculate improved shear correction factors

        Reference:

            Vlachoutsis, S. "Shear correction factors for plates and shells",
            Int. Journal for Numerical Methods in Engineering, Vol. 33,
            1537-1552, 1992.

            http://onlinelibrary.wiley.com/doi/10.1002/nme.1620330712/full


        Using "one shear correction factor" (see reference), assuming:

        - constant G13, G23, E1, E2, nu12, nu21 within each ply
        - g1 calculated using z at the middle of each ply
        - zn1 = Laminate.offset

        Returns
        -------

        k13, k23 : tuple
            Shear correction factors. Also updates attributes: `scf_k13`
            and `scf_k23`.