anchor.diff:          rpar(5) = kfn
           rpar(6) = R_a
   
  -        call dop853(n,structure_factor_phonon,kn_start,y,kn_end,h, &
  +        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)
  
positive.code:         itol = 0
          rtol = 1.0e-4
          atol = 1.0e-5
          iout = 0    ! want dense output
          lout = error_unit
          lipar = num_ipar
          lrpar = num_rpar
          n = num_integration_variables
                    
          y(1) = 0.0
  
  
negative.code:         real(dp), intent(in) :: nn, nion, temperature
          type(composition_info_type), intent(in) :: ionic
          integer, intent(out) :: ierr
          real(dp) :: nu      ! neutron-impurity scattering frequency
  
          integer, parameter :: num_integration_variables = 1
          integer, parameter :: num_ipar = 0
          integer, parameter :: num_rpar = 6
      
          real(dp) :: ne      ! number density of electrons
          real(dp) :: kFn     ! neutron Fermi wavevector 
          real(dp) :: EFn     ! neutron Fermi energy
          real(dp) :: V       ! effective neutron-impurity potential 
          real(dp) :: Lambda_n_imp    ! Coulomb logarithm (Potekhin et al. 1999)
          real(dp) :: fac, sf_frac, Tc, xFn
  