try_to_optimise_passing_grid Subroutine

  subroutine try_to_optimise_passing_grid(optimal_split, optimal_points, &
       use_max_error, noptimal_points)
    use gauss_quad, only: get_legendre_grids_from_cheb
    use theta_grid, only: ntgrid, bmag, bmax
    use splines, only: new_spline, spline, splint
    use optionals, only: get_option_with_default
    implicit none
    real, intent(out) :: optimal_split
    integer, intent(out) :: optimal_points
    logical, intent(in), optional :: use_max_error
    integer, intent(in), optional :: noptimal_points
    real, dimension(:), allocatable :: lambda
    real, dimension(:, :), allocatable :: weights
    real, dimension (:), allocatable :: wx_tmp, xx_tmp
    real, dimension(:), allocatable :: analytic_sum, b_ratio
    real, dimension(:), allocatable :: error
    real, dimension(:), allocatable :: error_history
    real, dimension(*), parameter :: trial_splits = &
         [0.001, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1, 0.5]
    integer, parameter :: ntrial = size(trial_splits)
    real, parameter :: min_split = trial_splits(1), max_split = trial_splits(ntrial)
    integer :: high_ntrial, itrial, il
    type(spline) :: the_spline
    integer :: minimum_location
    logical :: use_max_error_local
    logical, parameter :: debug = .false.
    integer, parameter :: high_res_factor = 100
    real :: error_from_no_split
    use_max_error_local = get_option_with_default(use_max_error, .false.)

    allocate (analytic_sum(-ntgrid:ntgrid))
    allocate (b_ratio(-ntgrid:ntgrid))
    allocate (error(-ntgrid:ntgrid))

    b_ratio = bmax / bmag
    analytic_sum = 2.0 * (1 - sqrt(b_ratio - 1)/sqrt(b_ratio))

    ! Determine how many points to use in the uppper region.
    ! Here we set it to about half if not passed in.
    ! For now this is a user choice/hard-coded value. We could
    ! imagine also optimising this value.
    optimal_points = get_option_with_default(noptimal_points, ng2 - (1 + ng2/2))
    if (optimal_points <= 0) optimal_points = ng2/2

    allocate(error_history(ntrial))
    allocate(xx_tmp(ng2), wx_tmp(ng2), lambda(ng2))
    allocate(weights(-ntgrid:ntgrid, ng2))

    ! Determine the initial error metric if we don't split
    call get_legendre_grids_from_cheb (1., 0., xx_tmp, wx_tmp)
    lambda = (1.0 - xx_tmp**2)/bmax
    ! Transform the Legendre weights to the lambda grid
    do il = 1, ng2
       weights(:,il) = wx_tmp(il)*2.0*sqrt((bmag/bmax) &
            *((1.0/bmax-lambda(il))/(1.0/bmag-lambda(il))))
    end do

    error = analytic_sum - sum(weights, dim=2)

    error_from_no_split = get_error_metric(error)

    ! Run with a few splits
    do itrial = 1, ntrial
       call get_legendre_grids_from_cheb (1., trial_splits(itrial), xx_tmp(:optimal_points), wx_tmp(:optimal_points))
       call get_legendre_grids_from_cheb (trial_splits(itrial), 0.0, xx_tmp(1+optimal_points:), wx_tmp(1+optimal_points:))

       lambda = (1.0 - xx_tmp**2)/bmax

       ! Transform the Legendre weights to the lambda grid
       do il = 1, ng2
          weights(:,il) = wx_tmp(il)*2.0*sqrt((bmag/bmax) &
               *((1.0/bmax-lambda(il))/(1.0/bmag-lambda(il))))
       end do

       error = analytic_sum - sum(weights, dim=2)

       error_history(itrial) = get_error_metric(error)
    end do

    deallocate(error)
    deallocate(xx_tmp)

    if (debug) print*,"History",error_history

    ! Spline the error metric vs splits
    the_spline = new_spline(trial_splits, error_history)

    high_ntrial = ntrial * high_res_factor
    allocate(error(high_ntrial))
    allocate(xx_tmp(high_ntrial))

    ! Evaluate the spline on a high resolution grid
    do itrial = 1, high_ntrial
       xx_tmp(itrial) = min_split + (itrial-1)*(max_split-min_split)/(high_ntrial-1)
       error(itrial) = splint(xx_tmp(itrial), the_spline)
    end do

    ! Choose the location of the minimum error metric as the optimal split point
    minimum_location = minloc(abs(error), dim = 1)

    ! We should probably calculate the actual error with this split
    ! choice in case our spline is not a good representation of the
    ! data. If the real error is larger than any of the real measurements
    ! we could just opt to use that instead?
    optimal_split = xx_tmp(minimum_location)

    if (debug) print*,"Optimal split : ",optimal_split," with error ",error(minimum_location)
    if (minval(abs(error)) > error_from_no_split) then
       optimal_split = -1
       optimal_points = -1
       if (debug) print*,'Minimum error from split larger than no split'
    end if
  contains
    pure real function get_error_metric(errors) result(metric)
      implicit none
      real, dimension(:), intent(in) :: errors
       ! Choose error metric we will later minimise
       if (use_max_error_local) then
          ! Maximum error
          metric = maxval(abs(errors))
       else
          ! Average error
          metric = sum(abs(errors))/size(errors)
       end if
    end function get_error_metric

  end subroutine try_to_optimise_passing_grid