anchor.diff:  ! Input arguments
       integer nt, np ! number of observations, planets
       real (kind=8), dimension(nt) :: time, vel
  -    real (kind=8), dimension(np) :: p, k, ecc, omega, t0
  +    real (kind=8), dimension(np) :: P, K, ecc, omega, t0
       real (kind=8) :: vsys
   
   ! Local variables
  -    !integer,parameter :: sp = selected_real_kind(p=6,r=37)
  -    integer,parameter :: dp = selected_real_kind(p=15,r=307)
       integer :: i
  -
  -    real(dp), parameter :: pi = 3.1415926535897932384626433832795029_dp
  -    real(dp), parameter :: twopi = 2.0_dp * pi 
  +    real(kind=8), parameter :: pi = 3.1415926535897932384626433832795029d0
  +    real(kind=8), parameter :: twopi = 2.0d0 * pi 
   
   !     vel = vsys
       do i=1,np
  
positive.code:   end function rv_curve
  
  
    elemental real(dp) function true_anom(M, ecc) result(nu)
      real(dp), intent(in) :: M, ecc
      real(dp) :: E, cosE
  
      E = ecc_anom(M, ecc)
  
  
negative.code: 
  contains
    
    elemental real(kind=8) function rv_curve(time, P, K, ecc, omega, t0) result(vel)
    ! RV curve for one planet - see Balan & Lahav (2009), MNRAS, 394, 1936
      real(kind=8), intent(in) :: time
      real(kind=8), intent(in) :: P, K, ecc, omega, t0
      real(kind=8) :: M, f 
  
      M = twopi * (time-t0)/P
      f = true_anom(M, ecc)
  