anchor.diff:    !****************************************************
     ! linear combination of stiffness and mass matrix
     !****************************************************
  -  function S_S(alpha, a, b)
  +  function S_S(alpha)
   
       type(body)   :: alpha
  -    real(dp)     :: a, b
       type(matrix) :: S_S
   
  -    S_S = alpha%K  + b*alpha%M ! K + b M
  +    S_S = alpha%K  + bb*alpha%M ! K + b M
   
     end function S_S
   
     !****************************************************
     !****************************************************
  -  function S_RS(alpha, a, b)
  +  function S_RS(alpha)
   
       type(body)   :: alpha
  -    real(dp)     :: a, b
       type(matrix) :: S_RS(1,4)
       type(matrix) :: O ! zero matrix
   
  
positive.code: 
    end function jac_full
  
    !****************************************************
    ! returns the rigid body terms in the jacobian matrix
    !****************************************************
    function jac_rigid(alpha, a)
  
      type(body)      :: alpha
      real(dp)        :: a
      real(dp)        :: jac_rigid(12,12)
      type(matrix)    :: DR1 (4,4)
  
      ! assemble the jacobian of type(matrix)
      DR1              = D_R(alpha, a)   
  
      ! now return in real-matrix form needed for computations
      jac_rigid       = matrix(DR1)
  
    end function jac_rigid
  
    !*************************************************************************
    ! find the jacobian of the equation of motion for the supplied BODY alpha
    !*************************************************************************
    function D_R(alpha, a)
      use dispmodule
      type(body)    :: alpha
      real(dp)      :: a
      type(matrix)  :: D_R(4,4) ! where 4 is the state vector length
      type(matrix)  :: O ! zero matrix
      !    type(matrix)  :: U ! unit matrix
  
negative.code: 
      S_RS(1,1) = O
      S_RS(1,2) = O
      S_RS(1,3) = bb*alpha%p + aa*skew(alpha%omega)*alpha%p !bp +a wx p
      S_RS(1,4) = bb*alpha%h + aa*(skew(alpha%v)*alpha%p  + skew(alpha%omega)*alpha%h) !bh +a(vx p + wx h)
  
    end function S_RS
  
    !****************************************************
    !****************************************************
    function S_SR(alpha)
  
      type(body)   :: alpha
      type(matrix) :: S_SR(4,1)
      type(matrix) :: O ! zero matrix
  
  