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:     I = matrix(eye(num_spat_dim))
  
      ! assemble jacobian
      D_R(1,1) = a * alpha%C_mat
      D_R(2,1) = O
      D_R(3,1) = O
      D_R(4,1) = O
  
      D_R(1,2) = skew(alpha%C_mat * alpha%r_dot) * alpha%S
      D_R(2,2) = alpha%S_dot + skew(alpha%S*alpha%theta_dot)*alpha%S + a * alpha%S
      D_R(3,2) = O
      D_R(4,2) = O
  
      D_R(1,3) = -1.0_dp*I
      D_R(2,3) = O
      D_R(3,3) = alpha%mass*(a*I + skew(alpha%omega))
      D_R(4,3) = a*skew(alpha%c) + (skew(alpha%c) * skew(alpha%omega))
  
      D_R(1,4) = O
      D_R(2,4) = -1.0_dp*I
      D_R(3,4) = -a*skew(alpha%c) + skew(skew(alpha%c)*alpha%omega)  - alpha%mass*skew(alpha%V) - skew(alpha%omega)*skew(alpha%c)
      D_R(4,4) = a*alpha%J  + skew(alpha%omega)*alpha%J - skew(alpha%J*alpha%omega) -skew(alpha%c)*skew(alpha%V)
  
    end function D_R
  
    !**********************************************
    ! jacobian of the structural degree of freedom
    !**********************************************
    function S_R(alpha, a)
  
      type(body)   :: alpha
      real(dp)     :: a
      type(matrix) :: S_R(4,4) ! where 4 is the state vector length
      type(matrix)  :: O ! zero matrix
  
  
negative.code: 
      S_SR(1,1) = O
      S_SR(2,1) = O
      S_SR(3,1) = aa * trans(alpha%p) !a p^T
      S_SR(4,1) = aa * trans(alpha%p) !a h^T
  
    end function S_SR
  
  