anchor.diff:      	nextra = 0
       	if (nest_context / 100 == 2) then
       		using_gp = .true.
  -    		nextra = 0  ! hyperparameters of GP
  +    		nextra = 1  ! hyperparameters of GP
       	end if
       end if
  -
  -    !! some parameters depend on the ones set on the namelist
  -    !nplanets = sdim / 5 ! this is too restrictive to the extra parameters
  -    nplanets = mod(nest_context,100) / 10 ! this is good for now
  -    allocate(spriorran(sdim, 2))
  -    allocate(nest_pWrap(sdim))
  -    nest_nPar=sdim
   
   	!doing debug?
   	doing_debug = .false.
   
   
  -	!load data
  +	!! load data
   	open(unit=15, file='input.rv', status="old")
   	read(15, *) line ! header
   
  -	read(15, '(a)') line ! number of observatories
  +	read(15, '(a)') line
   	line = trim(line(index(line, "#")+1:))
  -	read(line, *) nobserv
  -
  -	allocate(n_each_observ(nobserv))
  -	read(15, '(a)') line ! number of measurements in each observatory
  +	read(line, *) nobserv ! number of observatories
  +    if (nobserv < 1) stop 'Error parsing input file'
  +
  +	allocate(n_each_observ(nobserv)) ! number of measurements in each observatory
  +	read(15, '(a)') line
   	line = trim(line(index(line, "#")+1:)) ! remove the #
   	delim = ' '
   	call parse(line, delim, args, nargs) ! split whitespaces
  
positive.code: 		spriorran(i:,1)= -kmax
  		spriorran(i:,2)= kmax
  
  	else
  		!! systematic velocity(ies), Uniform, -kmax - kmax
  		i = sdim-nextra+1 ! index of jitter parameter
  		spriorran(i:,1)= -kmax
  		spriorran(i:,2)= kmax
  
  	end if	
  
   	!stop "we have to stop here"
     	call nest_Sample
  
     	! deallocate memory
  
negative.code:     end if
  
      ! extra parameters are systematic velocities for each observatory 
      ! plus, if present, the jitter or the hyperparameters
      nextra = nextra + nobserv
  
  
      if (using_gp) then
          !write(*,*) 'Using Gaussian Process model'
      	!k3 = DiagonalKernel((/1.d0/))
      	!kernel_to_pass => k3
      	!k5 = ExpSquaredKernel((/ 1.d0, 1.d0 /))
      	!kernel_to_pass => k5
          k6 = ExpSineSquaredKernel((/ 1.d0, 25.d0 /))
          kernel_to_pass => k6
  
      	gp1 = GP(k6%evaluate_kernel(times, times), kernel_to_pass)
      	gp1%mean_fun => mean_fun_keplerian
  		gp_n_planets = 1
  		gp_n_planet_pars = 6
  