anchor.diff:  
       ! Deal with looping boundaries properly.
       if (EA_new .gt. twopi) then
  -      EA_new = 0.5 * (EA + twopi)
  -    elseif (EA_new .lt. 0.0) then
  -      EA_new = 0.5 * EA 
  +      EA_new = 0.5d0 * (EA + twopi)
  +    elseif (EA_new .lt. 0.0d0) then
  +      EA_new = 0.5d0 * EA 
       endif
   
     end function halley
           
  -  elemental real(dp) function strict_iteration(EA, e, M) result(EA_new)
  +  elemental real(kind=8) function strict_iteration(EA, e, M) result(EA_new)
     ! Solve kepler's equation by a simple iteration scheme
  -    real(dp), intent(in) :: EA, e, M
  +    real(kind=8), intent(in) :: EA, e, M
       EA_new = M+e*sin(EA)
     end function strict_iteration
           
  -  elemental real(dp) function newton(EA, e, M) result(EA_new)
  +  elemental real(kind=8) function newton(EA, e, M) result(EA_new)
     ! The following routine is used to solve kepler's equation using
     ! Newton's method. It is very fast and reliable for small values of
     ! e, but can be wildly erratic for e close to 1. See, e.g, the discussion
     ! in Jean Meeus, Astronomical Algorithms, Willmann-Bell, 1991, 181-193.
  -    real(dp), intent(in) :: EA, e, M
  -    EA_new = EA + (M + e*sin(EA) - EA)/(1 - e*cos(EA))
  +    real(kind=8), intent(in) :: EA, e, M
  +    EA_new = EA + (M + e*sin(EA) - EA)/(1d0 - e*cos(EA))
     end function newton
     
     
  
positive.code:     endif
  
      iter = 0
      EA_old = 0.5_dp * pi 
      if (ecc < 0.8d0) then
        EA_new = newton(EA_old, ecc, M)
        do while (abs(EA_old - EA_new) >= derror .and. iter<200)
  
negative.code:   end function ecc_anom
          
          
    real(kind=8) function halley(EA, e, M) result(EA_new)
      real(kind=8), intent(in) :: EA, e, M
      real(kind=8) :: spsi0, f, fp, fpp
  
      ! Compute the function and derivatives.
      spsi0 = sin(EA)
      f = EA - e * spsi0 - M
      fp = 1.d0 - e * cos(EA)
      fpp = e * spsi0
  
      ! Take a second order step.
  