anchor.diff:  
     function find_collision_time(a_i, a_j, R, t0, t1) BIND(C)
       implicit none
  -    double precision :: find_collision_time
  -    double precision, dimension(3,3), intent(in) :: a_i, a_j
  -    double precision, value, intent(in) :: R, t0, t1
  -    double precision, dimension(3,3) :: a_ji
  -    double precision, dimension(5) :: P
  -    double precision :: a_x, a_y, b_x, b_y, c_x, c_y
  +    real(c_double) :: find_collision_time
  +    real(c_double), dimension(3,3), intent(in) :: a_i, a_j
  +    real(c_double), value, intent(in) :: R, t0, t1
  +    real(c_double), dimension(3,3) :: a_ji
  +    real(c_double), dimension(5) :: P
  +    real(c_double) :: a_x, a_y, b_x, b_y, c_x, c_y
       a_ji = a_i - a_j
       a_x = a_ji(1,3)
       a_y = a_ji(3,3)
  
positive.code: 
    function sort_complex_conjugate_pairs(roots) BIND(C)
      implicit none
      integer :: sort_complex_conjugate_pairs
      double complex, dimension(4), intent(inout) :: roots
      integer :: npairs = 0
      integer :: i, j
      double complex :: r, r_conj
      i = 1
      do while (i <= 3)
         r = roots(i)
         if (abs(DIMAG(r)) > IMAG_TOLERANCE) then
            do j = 1, 3-i
               r_conj = roots(i+j)
               if (       abs(DREAL(r) - DREAL(r_conj)) < ZERO_TOLERANCE &
                    .and. abs(DIMAG(r) + DIMAG(r_conj)) < ZERO_TOLERANCE) then
  
negative.code:        endif
         S = 0.5*sqrt(SSx4)
      endif
      sqrtp = sqrt(-SSx4 - 2*p + q/S)
      sqrtm = sqrt(-SSx4 - 2*p - q/S)
      out(1) = -0.25*b - S + 0.5*sqrtp;
      out(2) = -0.25*b - S - 0.5*sqrtp;
      out(3) = -0.25*b + S + 0.5*sqrtm;
  