anchor.diff:      O = matrix(zeros(num_spat_dim))
   
       S_R(1,1) = O
  +    S_R(1,2) = O
  +    S_R(1,3) = O
  +    S_R(1,4) = O
  +
       S_R(2,1) = O
  +    S_R(2,2) = O
  +    S_R(2,3) = O
  +    S_R(2,4) = O
  +
       S_R(3,1) = O
  +    S_R(3,2) = O
  +    S_R(3,3) = O
  +    S_R(3,4) = - skew(alpha%p*alpha%qs_dot)
  +
       S_R(4,1) = O
  -
  -    S_R(1,2) = O
  -    S_R(2,2) = O
  -    S_R(3,2) = O
       S_R(4,2) = O
  -
  -    S_R(1,3) = O
  -    S_R(2,3) = O
  -    S_R(3,3) = O
       S_R(4,3) = - skew(alpha%p*alpha%qs_dot)
  -
  -    S_R(1,4) = O
  -    S_R(2,4) = O
  -    S_R(3,4) = - skew(alpha%p*alpha%qs_dot)
       S_R(4,4) = - skew(alpha%h*alpha%qs_dot)
   
     end function S_R
  
positive.code: 
      type(body)   :: alpha
      real(dp)     :: a
      type(matrix) :: S_SR(1,4)
      type(matrix) :: O ! zero matrix
  
      O = matrix(zeros(num_spat_dim))
  
      S_SR(1,1) = O
      S_SR(1,2) = O
      S_SR(1,3) = a * trans(alpha%p) !a p^T
      S_SR(1,4) = a * trans(alpha%p) !a h^T
  
    end function S_SR
  
  
negative.code: 
      type(body)   :: alpha
      real(dp)     :: a, b
      type(matrix) :: S_RS(1,4)
      type(matrix) :: O ! zero matrix
  
      O = matrix(zeros(num_spat_dim))
  
      S_RS(1,1) = O
      S_RS(1,2) = O
      S_RS(1,3) = b*alpha%p + a*skew(alpha%omega)*alpha%p !bp +a wx p
      S_RS(1,4) = b*alpha%h + a*(skew(alpha%v)*alpha%p  + skew(alpha%omega)*alpha%h) !bh +a(vx p + wx h)
  
    end function S_RS
  
  