anchor.diff:          itol = 0
           rtol = 1.0e-4
           atol = 1.0e-5
  -        iout = 0    ! want dense output
  +        iout = 0    ! solout is never called for output
           lout = error_unit
           lipar = num_ipar
           lrpar = num_rpar
  -        n = num_integration_variables
                       
           !calculate n_in and R_a using Sly4 model for now
           isospin = 1.0 - 2.0*ionic%Z/ionic%A
  
positive.code:         real(dp) :: qD, qi, ktf, beta, qs, wf, sf, L2, eta_n_0, G_k, Tp, eta_n
          real(dp) :: isospin, n_0, n_2, n_l, n_in
          real(dp), dimension(:), pointer :: y
          integer, intent(out) :: ierr
          integer ,dimension(:), pointer :: iwork => null()
          real(dp), dimension(:), pointer :: work => null()
          real(dp), dimension(1) :: rtol, atol
          integer, dimension(:), pointer :: ipar => null()
          real(dp), dimension(:), pointer :: rpar => null()
          real(dp) :: h, kn_end, kn_start, kn, kfi
          integer :: liwork, lwork, itol, lipar, lrpar
          integer :: n, idid, lout, iout, npts, num_integration_variables, num_ipar, num_rpar
  
          allocate(y(1))
          num_ipar = 0
          num_rpar = 6    
          kFn = (threepisquare*nn)**onethird
          kfi = (threepisquare*nn/ionic%Yn*(1.0-ionic%Yn)/ionic%A)**onethird
          num_integration_variables = 1
          call dop853_work_sizes( &
          &   num_integration_variables,num_integration_variables,liwork,lwork)
          allocate(iwork(liwork), work(lwork),ipar(num_ipar), rpar(num_rpar))
  
negative.code:         rpar(5) = kfn
          rpar(6) = R_a
  
          call dop853(num_integration_variables, &
          &   structure_factor_phonon,kn_start,y,kn_end,h, &
          &   int_default_max_step_size,int_default_max_steps, &
          &   rtol,atol,itol, null_solout, iout, work, lwork, iwork, liwork,  &
          &   num_rpar, rpar, num_ipar, ipar, lout, idid)
  