diff --git a/moorpy/line.py b/moorpy/line.py
index 857defb..b0e7b58 100644
--- a/moorpy/line.py
+++ b/moorpy/line.py
@@ -851,7 +851,7 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0):
         self.TB = np.linalg.norm(self.fB)
         
         # Compute transverse (out-of-plane) stiffness term
-        if fAV > fAH:  # if line is more vertical than horizontal, 
+        if LH < 0.01*abs(LV):  # if line is nearly vertical (note: this theshold is unverified)
             Kt = 0.5*(fAV-fBV)/LV  # compute Kt based on vertical tension/span
         else:  # otherwise use the classic horizontal approach
             Kt = -fBH/LH
@@ -860,7 +860,7 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0):
         # save 3d stiffness matrix in global orientation for both line ends (3 DOF + 3 DOF)
         self.KA  = from2Dto3Drotated(info['stiffnessA'],  Kt, R.T)  # reaction at A due to motion of A
         self.KB  = from2Dto3Drotated(info['stiffnessB'],  Kt, R.T)  # reaction at B due to motion of B
-        self.KBA = from2Dto3Drotated(info['stiffnessBA'], Kt, R.T)  # reaction at B due to motion of A
+        self.KBA = from2Dto3Drotated(info['stiffnessBA'],-Kt, R.T)  # reaction at B due to motion of A
         
         
         # ----- calculate current loads if applicable, for use next time -----
diff --git a/moorpy/subsystem.py b/moorpy/subsystem.py
index b529cbe..f251a7f 100644
--- a/moorpy/subsystem.py
+++ b/moorpy/subsystem.py
@@ -249,6 +249,7 @@ def staticSolve(self, reset=False, tol=0.01, profiles=0):
         # outputs should be pointList[0] and [N] .r
         dr =  self.rB - self.rA
         LH = np.hypot(dr[0], dr[1])         # horizontal spacing of line ends
+        LV = dr[2]                          # vertical rise from end A to B
         self.pointList[ 0].setPosition([ -self.span   , 0, self.rA[2]])
         self.pointList[-1].setPosition([ -self.span+LH, 0, self.rB[2]])
         
@@ -273,17 +274,18 @@ def staticSolve(self, reset=False, tol=0.01, profiles=0):
         # save end forces and stiffness matrices (first in local frame)
         self.fA_L = self.pointList[ 0].getForces(xyz=True) # force at end A
         self.fB_L = self.pointList[-1].getForces(xyz=True) # force at end B
-        '''
-        self.KA_L  = K[:3,:3]                # reaction at A due to motion of A
-        self.KB_L  = K[3:,3:]                # reaction at B due to motion of B
-        self.KAB_L = K[3:,:3]                # reaction at B due to motion of A
-        '''
+        
+        # Compute transverse (out-of-plane) stiffness term
+        if LH < 0.01*abs(LV):  # if line is nearly vertical (note: this theshold is unverified)
+            Kt = 0.5*(self.fA_L[2] - self.fB_L[2])/LV  # compute Kt based on vertical tension/span
+        else:  # otherwise use the classic horizontal approach
+            Kt = -self.fB_L[0]/LH  
         
         # expand to get 3D stiffness matrices
         R = np.eye(3)
-        self.KA_L  = from2Dto3Drotated(K[:2,:2], -self.fB_L[0], LH, R.T)  # reaction at A due to motion of A
-        self.KB_L  = from2Dto3Drotated(K[2:,2:], -self.fB_L[0], LH, R.T)  # reaction at B due to motion of B
-        self.KBA_L = from2Dto3Drotated(K[2:,:2], -self.fB_L[0], LH, R.T)  # reaction at B due to motion of A
+        self.KA_L  = from2Dto3Drotated(K[:2,:2],  Kt, R.T)  # reaction at A due to motion of A
+        self.KB_L  = from2Dto3Drotated(K[2:,2:],  Kt, R.T)  # reaction at B due to motion of B
+        self.KBA_L = from2Dto3Drotated(K[2:,:2], -Kt, R.T)  # reaction at B due to motion of A
         
         self.TA = np.linalg.norm(self.fA_L)  # tensions [N]
         self.TB = np.linalg.norm(self.fB_L)
diff --git a/tests/test_system.py b/tests/test_system.py
index f071e4c..81d3882 100644
--- a/tests/test_system.py
+++ b/tests/test_system.py
@@ -340,6 +340,44 @@ def test_seabed(CB):
                     np.hstack([ms2.pointList[0].getForces(), ms2.pointList[-1].getForces()]), rtol=0, atol=10.0, verbose=True)
 
 
+
+def test_6DOF_stiffness():
+    '''Tests 6x6 stiffness matrix of a complex system to check off diagonals.'''
+    
+    # Create new MoorPy System and set its depth
+    ms = mp.System(file='tests/case8.dat')
+
+    # ----- run the model to demonstrate -----
+
+    ms.initialize()                                             # make sure everything's connected
+
+    # find equilibrium of just the mooring lines (the body is 'coupled' and will by default not be moved)
+    ms.solveEquilibrium(tol=0.0001)                                       # equilibrate
+
+
+    #Kmp  = ms.getSystemStiffness( DOFtype='both', dx=0.02, dth=0.001)
+    #KmpA = ms.getSystemStiffnessA(DOFtype='both', )
+    
+    Kn = ms.getCoupledStiffness()
+    Ka = ms.getCoupledStiffnessA()
+    
+
+    #Kn = Kmp
+    #Ka = KmpA
+    #Kn = ms.bodyList[0].getStiffness(tol=0.00001, dx = 0.001)
+    #Ka = ms.bodyList[0].getStiffnessA(lines_only=True)
+    
+    print(Kn)
+    print(Ka)
+    print(Ka-Kn)
+    #print(Ka/Kn)
+    
+    assert_allclose(Ka[:3,:3], Kn[:3,:3], rtol=0.02, atol=5e3, verbose=True) # translational
+    assert_allclose(Ka[3:,:3], Kn[3:,:3], rtol=0.02, atol=5e4, verbose=True) #
+    assert_allclose(Ka[:3,3:], Kn[:3,3:], rtol=0.02, atol=5e4, verbose=True) #
+    assert_allclose(Ka[3:,3:], Kn[3:,3:], rtol=0.02, atol=2e6, verbose=True) # rotational
+
+
 if __name__ == '__main__':
     
     #test_basic()