anchor.diff:    end function rv_curve
   
   
  -  elemental real(dp) function true_anom(M, ecc) result(nu)
  -    real(dp), intent(in) :: M, ecc
  -    real(dp) :: E, cosE
  +  elemental real(kind=8) function true_anom(M, ecc) result(nu)
  +    real(kind=8), intent(in) :: M, ecc
  +    real(kind=8) :: E, cosE
   
       E = ecc_anom(M, ecc)
   
  
positive.code: 
  contains
    
    elemental real(dp) function rv_curve(time, P, K, ecc, omega, t0) result(vel)
    ! RV curve for one planet - see Balan & Lahav (2009), MNRAS, 394, 1936
      real(dp), intent(in) :: time
      real(dp), intent(in) :: P, K, ecc, omega, t0
      real(dp) :: M, f 
  
      M = twopi * (time-t0)/P
      f = true_anom(M, ecc)
  
negative.code:   
  
  
    elemental real(kind=8) function ecc_anom(M, ecc) result(E)
      real(kind=8), intent(in) :: M, ecc
      real(kind=8), parameter :: derror = 1.0d-7
      integer :: iter
      real(kind=8) :: EA_new, EA_old
      
      ! catch special case
      if (ecc .eq. 0.d0) then
  