module field_checks
    use constants, only: dp, pi
    use field_base, only: field_t

    implicit none
    private

    public :: suspect_omnigenous_origin_not_minimum

contains

    function suspect_omnigenous_origin_not_minimum(field, N_tor, M_pol, retol)
        use find_extrema, only: find_local_minima
        use find_extrema, only: find_global_extrema
        use field_base, only: field_t
        use fieldline_mod, only: fieldline_t
        use field_along_fieldline, only: set_field_and_fieldline
        use field_along_fieldline, only: unset_field_and_fieldline
        use field_along_fieldline, only: B_mod_along_fieldline, dB_dphi_along_fieldline
        class(field_t), intent(in) :: field
        real(dp), intent(in) :: N_tor, M_pol
        real(dp), intent(in), optional :: retol
        logical :: suspect_omnigenous_origin_not_minimum

        real(dp) :: iota
        real(dp), dimension(2) :: chi_interval, phi_interval
        type(fieldline_t) :: fieldline
        real(dp), dimension(:), allocatable :: phi, phi_error, B, dB_dphi
        integer :: n_min, n_extrema, idx
        real(dp) :: phi_max, phi_tol
        real(dp), parameter :: height_retol = 0.3_dp
        real(dp), dimension(2) :: extrema
        real(dp) :: B_min, B_max, B_range, B_at_origin, height
        real(dp) :: B_min_error, B_max_error, B_range_error
        real(dp) :: height_error

        ! as M_pol and N_tor are whole numbers that must no both be zero
        ! M_pol*pi - N_tor should never zero
        if (abs(M_pol*pi - N_tor) < 1e-8) then
            print *, "Error: (M_pol*pi - N_tor) must not be (close) zero."
            print *, "abs(M_pol*iota - N_tor) = ", abs(M_pol*pi - N_tor)
            error stop
        end if

        iota = pi
        fieldline%iota = iota
        fieldline%phi_0 = 0.0_dp
        fieldline%theta_0 = 0.0_dp

        suspect_omnigenous_origin_not_minimum = .false.

        chi_interval = [-pi, pi]
        phi_interval = chi_interval/(M_pol*iota - N_tor)

        call set_field_and_fieldline(field, fieldline)

        call find_local_minima(B_mod_along_fieldline, phi_interval, phi, phi_error)
        n_min = size(phi)
        if (n_min == 0) then
            print *, "Warning in suspect_omnigenous_origin_not_minimum: ", &
                "No local dB_dphifound in well around origin! "
            suspect_omnigenous_origin_not_minimum = .true.
            call unset_field_and_fieldline()
            return
        end if
        allocate (B(n_min))
        allocate (dB_dphi(n_min))
        call B_mod_along_fieldline(phi, B)
        call dB_dphi_along_fieldline(phi, dB_dphi)
        idx = minloc(B, dim=1)
        B_min_error = abs(dB_dphi(idx))*phi_error(idx)
        B_min = B(idx)
        deallocate (B, dB_dphi)

        phi_tol = abs(phi_interval(2) - phi_interval(1))*1e-3_dp
        extrema = find_global_extrema(B_mod_along_fieldline, phi_interval, &
                                      abstol=phi_tol)
        n_extrema = size(extrema)
        allocate (B(n_extrema))
        allocate (dB_dphi(n_extrema))
        call B_mod_along_fieldline(extrema, B)
        call dB_dphi_along_fieldline(extrema, dB_dphi)
        idx = maxloc(B, dim=1)
        B_max_error = abs(dB_dphi(idx))*phi_tol
        B_max = B(idx)
        deallocate (B, dB_dphi)

        call unset_field_and_fieldline()

        ! is the B-difference of origin and minimum significant
        ! compared to the B-range?
        B_range = B_max - B_min
        B_range_error = B_min_error + B_max_error
        call field%compute_B_mod(0.0_dp, 0.0_dp, B_at_origin)
        height = B_at_origin - B_min
        height_error = B_min_error
        if (present(retol)) height_error = height_error + retol*B_at_origin
        if ((height - height_error) > height_retol*(B_range + B_range_error)) then
            print *, "Detected that B at origin of provided field is"
            print *, "significantly above the minimum B i.e. difference > "
            print *, height_retol*100.0_dp, "% of the total B range!"
            print *, "(B_at_origin - B_min) = ", height
            print *, "with estimated error = ", height_error
            print *, "B_max = ", B_max, " B_min = ", B_min
            print *, "(B_max-B_min) = ", B_range
            print *, "with estimated error = ", B_range_error
            suspect_omnigenous_origin_not_minimum = .true.
        end if

    end function suspect_omnigenous_origin_not_minimum

end module field_checks
