setup_passing_lambda_grids Subroutine

  subroutine setup_passing_lambda_grids(lambda_grid, weights)
    use gauss_quad, only: get_legendre_grids_from_cheb, get_radau_gauss_grids
    use theta_grid, only: ntgrid, bmag, bmax
    implicit none
    !> The lambda grid points. Note this is intent `in out` as we only
    !> set a portion of the grid so may want to keep the rest of the
    !> input values.
    real, dimension(:), intent(in out) :: lambda_grid
    !> The integration weights for the lambda grid. This is intent `in out`
    !> for the same reason as `lambda_grid`.
    real, dimension(-ntgrid:, :), intent(in out) :: weights
    real, dimension (:), allocatable :: wx, xx_radau
    integer :: ig, il, passing_split_index
    real :: passing_split_value

    ! Note this is a module level quantity that we keep around for possible
    ! future use. For iproc /= 0 this allocation is actually done in [[broadcast_results]]
    if (.not. allocated(xx)) allocate (xx(ng2))

    if(radau_gauss_grid .and. grid_has_trapped_particles() ) then

       ! This grid uses a fixed endpoint such that il = ng2+1 has a finite weight at bounce points
       ! which contributes to the integration over passing particles
       ! which is why xx, wx are one element longer than the standard Legendre Gauss grid
       ! This grid should only be used when nlambda > ng2, i.e. there is a wfb and trapped particles

       allocate (xx_radau(ng2 +1))
       allocate (wx(ng2 +1))
       call get_radau_gauss_grids(0.,1., xx_radau, wx,report_in=.false.)
       xx = xx_radau(:ng2)
    else
       allocate (wx(ng2))
       if (split_passing_region) then
          call try_to_optimise_passing_grid(passing_split_value, passing_split_index)
          if (passing_split_index > 0) then
             call get_legendre_grids_from_cheb (1., passing_split_value, &
                  xx(:passing_split_index), wx(:passing_split_index))
             call get_legendre_grids_from_cheb (passing_split_value, 0.0, &
                  xx(1 + passing_split_index:), wx(1 + passing_split_index:))
          else
             call get_legendre_grids_from_cheb (1., 0., xx, wx)
          end if
       else
          call get_legendre_grids_from_cheb (1., 0., xx, wx)
       end if

    end if

    ! Store the location of the passing lambda grid points
    lambda_grid(:ng2) = (1.0 - xx(:ng2)**2)/bmax

    ! Transform the Legendre/Radau weights to the lambda grid
    do il = 1, ng2
       do ig = -ntgrid, ntgrid
          weights(ig,il) = wx(il)*2.0*sqrt((bmag(ig)/bmax) &
               *((1.0/bmax-lambda_grid(il))/(1.0/bmag(ig)-lambda_grid(il))))
       end do
    end do

    ! Assign the weight for wfb from the radau-gauss grid
    if (radau_gauss_grid .and. grid_has_trapped_particles() ) then
       il=ng2+1
       do ig = -ntgrid, ntgrid
          if (bmag(ig) < bmax) then ! i.e. we are not the bounce point for wfb
             ! Note this will currently always be exactly 0
             ! as lambda_grid(ng2+1) == 1.0/bmax
             weights(ig,il) = wx(il)*2.0*sqrt((bmag(ig)/bmax) &
                  *((1.0/bmax-lambda_grid(il))/(1.0/bmag(ig)-lambda_grid(il))))
          else ! at the bounce point
             weights(ig,il) = wx(il)*2.0
          endif
       end do
    end if

    if(allocated(wx)) deallocate(wx)
    if(allocated(xx_radau)) deallocate(xx_radau)
  end subroutine setup_passing_lambda_grids