anchor.diff:  
   
       !! Gaussian process variables
  -!     type(WhiteKernel), target :: k2
       type(DiagonalKernel), target :: k3
  -    !type(ExpKernel), target :: k4
       type(ExpSquaredKernel), target :: k5
       type(ExpSineSquaredKernel), target :: k6
       type(SumKernels), target :: k7
  
positive.code:     class(Kernel), pointer :: kernel_to_pass, subk1, subk2
  
      !! local variables
      real(kind=8), dimension(nt, nt) :: cov
  
      !WhiteKernel(3e-7) + 10. * ExpSquaredKernel(10) * ExpSine2Kernel(1, 30)
      !write(*,*) 'Using Gaussian Process model'
      k3 = DiagonalKernel((/3d-7/))
      k5 = ExpSquaredKernel((/ sqrt(10.d0) /))
      !kernel_to_pass => k5
      k6 = ExpSineSquaredKernel((/ 30.d0, sqrt(2.d0) /))
  
      k8 = ProductKernels((/1.d0, 1.d0 /))
  
negative.code: !     call gp1%sub_kernel1%get_kernel_pars(1)
  !     call gp1%sub_kernel1%get_kernel_pars(2)
  
      ! fill pred with systematic velocities
  !     do i=1,nobserv
  !         where(observ == i) pred = planetpar(size(planetpar)-)
  !     end do
      vel = planetpar(nPar)
  
      ! get the radial velocity model with these parameters (in vel)
      if (gp_n_planets > 0) then
          call get_rvN(times, &
                       planetpar(1:nPar-1:5), & ! periods for all planets
                       planetpar(2:nPar-1:5), & ! K for all planets
                       planetpar(3:nPar-1:5), & ! ecc for all planets
                       planetpar(4:nPar-1:5), & ! omega for all planets
                       planetpar(5:nPar-1:5), & ! t0 for all planets
                       0.d0, & ! systematic velocity
                       vel, nt, gp_n_planets)
      end if
  
      ! residuals: what is left when the planets (or just vsys) is subtracted from the data
      r = rvs - vel
  
      ! calculate predictions at observed times
      call gp1%predict(times, r, (/0.d0/), times, pred, cov, yerr=errors)
      pred = pred + vel
  
  
      return
  
  