anchor.diff:      I = matrix(eye(num_spat_dim))
   
       ! assemble jacobian
  -    D_R(1,1) = a * alpha%C_mat
  +    D_R(1,1) = aa * 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(2,2) = alpha%S_dot + skew(alpha%S*alpha%theta_dot)*alpha%S + aa * 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(3,3) = alpha%mass*(aa*I + skew(alpha%omega))
  +    D_R(4,3) = aa*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)
  +    D_R(3,4) = -aa*skew(alpha%c) + skew(skew(alpha%c)*alpha%omega)  - alpha%mass*skew(alpha%v) - skew(alpha%omega)*skew(alpha%c) 
  +    D_R(4,4) = aa*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)
  +  function S_R(alpha)
   
       type(body)   :: alpha
  -    real(dp)     :: a
       type(matrix) :: S_R(4,4) ! where 4 is the state vector length
       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:   !****************************************************
    ! linear combination of stiffness and mass matrix
    !****************************************************
    function S_S(alpha)
  
      type(body)   :: alpha
      type(matrix) :: S_S
  
      S_S = alpha%K  + bb*alpha%M ! K + b M
  
    end function S_S
  
    !****************************************************
    !****************************************************
    function S_RS(alpha)
  
      type(body)   :: alpha
      type(matrix) :: S_RS(1,4)
      type(matrix) :: O ! zero matrix
  
  