move_beta Subroutine

   pure subroutine move_beta( &
      np, beta, lower, upper, ssf, stpb, neta, eta, interval)
   !! Ensure range of bounds is large enough for derivative checking.
   !! Move beta away from bounds so that derivatives can be calculated.

      use odrpack_kinds, only: zero, one, three, ten, hundred

      integer, intent(in) :: np
         !! Number of parameters `np`.
      real(wp), intent(inout) :: beta(np)
         !! Function parameters.
      real(wp), intent(in) :: lower(np)
         !! !! Lower bound on `beta`.
      real(wp), intent(in) :: upper(np)
         !! Upper bound on `beta`.
      real(wp), intent(in) :: ssf(np)
         !! Scale used for the `beta`s.
      real(wp), intent(in) :: stpb(np)
         !! Relative step used for computing finite difference derivatives with respect to `beta`.
      integer, intent(in) :: neta
         !! Number of good digits in the function results.
      real(wp), intent(in) :: eta
         !! Relative noise in the function results.
      integer, intent(out) :: interval(np)
         !! Specifies which difference methods and step sizes are supported by the current
         !! interval `upper-lower`.

      ! Local scalars
      integer :: k
      real(wp) :: h, h0, h1, hc, hc0, hc1, stpl, stpr, typj

      ! VARIABLE DEFINITIONS (ALPHABETICALLY)
      !  H:        Relative step size for forward differences.
      !  H0:       Initial relative step size for forward differences.
      !  H1:       Default relative step size for forward differences.
      !  HC:       Relative step size for center differences.
      !  HC0:      Initial relative step size for center differences.
      !  HC1:      Default relative step size for center differences.
      !  K:        Index variable for BETA.
      !  STPL:     Maximum step to the left of BETA (-) the derivative checker will use.
      !  STPR:     Maximum step to the right of BETA (+) the derivative checker will use.
      !  TYPJ:     The typical size of the J-th unkonwn BETA.

      interval = 111
      do k = 1, np
         h0 = hstep(0, neta, 1, k, stpb, 1)
         hc0 = h0
         h1 = sqrt(eta)
         hc1 = eta**(one/three)
         h = max(ten*h1, min(hundred*h0, one))
         hc = max(ten*hc1, min(hundred*hc0, one))
         if (beta(k) == zero) then
            if (ssf(1) < zero) then
               typj = one/abs(ssf(1))
            else
               typj = one/ssf(k)
            end if
         else
            typj = abs(beta(k))
         end if
         stpr = (h*typj*sign(one, beta(k)) + beta(k)) - beta(k)
         stpl = (hc*typj*sign(one, beta(k)) + beta(k)) - beta(k)
         ! Check outer interval
         if (lower(k) + 2*abs(stpl) > upper(k)) then
            if (interval(k) >= 100) then
               interval(k) = interval(k) - 100
            end if
         elseif (beta(k) + stpl > upper(k) .or. beta(k) - stpl > upper(k)) then
            beta(k) = upper(k) - abs(stpl)
         elseif (beta(k) + stpl < lower(k) .or. beta(k) - stpl < lower(k)) then
            beta(k) = lower(k) + abs(stpl)
         end if
         ! Check middle interval
         if (lower(k) + 2*abs(stpr) > upper(k)) then
            if (mod(interval(k), 100) >= 10) then
               interval(k) = interval(k) - 10
            end if
         elseif (beta(k) + stpr > upper(k) .or. beta(k) - stpr > upper(k)) then
            beta(k) = upper(k) - abs(stpr)
         elseif (beta(k) + stpr < lower(k) .or. beta(k) - stpr < lower(k)) then
            beta(k) = lower(k) + abs(stpr)
         end if
         ! Check inner interval
         if (lower(k) + abs(stpr) > upper(k)) then
            interval(k) = 0
         elseif (beta(k) + stpr > upper(k)) then
            beta(k) = upper(k) - stpr
         elseif (beta(k) + stpr < lower(k)) then
            beta(k) = lower(k) - stpr
         end if
      end do

   end subroutine move_beta