module math
       implicit none

       real, private :: epsilon = 0.001
       integer, private :: maxloop = 100

       interface set_epsilon
          module procedure m_set_epsilon
       end interface

       interface set_maxloop
          module procedure m_set_maxloop
       end interface

       interface derivative
          module procedure m_derivative
       end interface

       interface newton
          module procedure m_newton
       end interface

       interface brent
          module procedure m_brent
       end interface

       interface bisection
          module procedure m_bisection
       end interface

       interface rootfind
          module procedure m_bisection
       end interface

       interface integral
          module procedure m_integral
       end interface

       interface trapz
          module procedure m_trapz
       end interface

       interface trapz_poisson_error
          module procedure m_trapz_poisson_error
       end interface

       interface sd
          module procedure m_sd
       end interface

contains

       subroutine m_set_epsilon(eps)
              real, intent(in) :: eps
              epsilon = abs(eps)
       end subroutine m_set_epsilon

       subroutine m_set_maxloop(max)
              integer, intent(in) :: max
              maxloop = max
       end subroutine m_set_maxloop

       real function m_derivative(f, x, dx0)
              implicit none

              real, intent(in) :: x
              real, intent(in), optional :: dx0
              real :: dx, prev_m_derivative
              integer :: i

              interface
                 real function f(x)
                        implicit none
                        real, intent(in) :: x
                 end function f
              end interface

              if (present(dx0)) then
                 dx = dx0
              else
                 dx = 100. * abs(epsilon*f(x))
              end if

              m_derivative = (f(x+dx/2.)-f(x-dx/2.))/dx

              do i=1,maxloop
                 ! write(*,*) " i=",i,"x=",x,"; dx=",dx,", f(x)=",f(x),", dfdx=",m_derivative
                 dx = dx/2.

                 if(dx<1e-7) stop "math::m_derivative: dx too small."

                 prev_m_derivative = m_derivative
                 m_derivative = (f(x+dx/2.)-f(x-dx/2.))/dx

                 if(abs((m_derivative-prev_m_derivative)/m_derivative)<epsilon) return
              end do
              stop "math::m_derivative: maxloop exceeded."
       end function m_derivative

       real recursive function m_newton(f, x, x0, loopnum, dx0) result (root)
              real, intent(in) :: x
              real, intent(in), optional :: dx0, x0
              real :: dx, dfdx, newx
              integer, optional :: loopnum
              integer :: l

              interface
                 real function f(x)
                        implicit none
                        real, intent(in) :: x
                 end function f
              end interface

              if(present(x0)) then
                 if( abs(x-x0) < epsilon) then
                    root = x
                    return
                 end if
              end if

              if (present(loopnum)) then
                 l = loopnum + 1
              else
                 l = 1
              end if

              if (l > maxloop) stop "math::m_newton: maximum number of loops exceeded."

              if (present(dx0)) then
                 dx = dx0
              else
                 dx = 10.
              end if

              dfdx = derivative(f, x, dx)
              newx = x - f(x)/dfdx

10            format (1(a," = ",i2),5("; ",a," = ",e10.4))

              !	write(*,10) "l",l,"x",x,"f(x)",f(x),"dfdx",dfdx,"newx",newx,"f(newx)",f(newx)
              root = m_newton(f, newx, x, l, dx)
       end function m_newton

       ! brent algorithm for finding roots (broken)
       real function m_brent(f, a0, b0)
              real, intent(in) :: a0, b0
              real :: a, b
              real :: dx
              integer :: i=0

              interface
                 real function f(x)
                        implicit none
                        real, intent(in) :: x
                 end function f
              end interface

              a = a0
              b = b0

              if(f(a)*f(b).lt.0.)then
                 m_brent=b
                 dx=-f(m_brent)/m_derivative(f, m_brent)
                 do while(abs(dx/m_brent).gt.epsilon .and. i.lt.maxloop)
                    i=i+1
                    dx=-f(m_brent)/m_derivative(f, m_brent)
                    m_brent=m_brent+dx
                    if(m_brent.gt.b.or.m_brent.lt.a) m_brent=(a+b)/2.     

                    if(F(a)*F(m_brent).gt.0.)then
                       a=m_brent           
                    else
                       b=m_brent
                    end if
                    write(*,*) i, a, b, dx, m_brent, m_derivative(f, m_brent)
                 end do
                 if(i.ge.maxloop) stop "math::m_brent: maximum number of loops exceeded."
              else
                 stop "math::m_brent: root not in interval."
              endif
       end function m_brent

       ! finds root by bisection. always works, inefficient.
       real recursive function m_bisection(func, a, b, loopnum) result(root)
              real, intent(in) :: a, b
              integer, intent(in), optional :: loopnum
              integer :: l

              interface
                 real function func(x)
                        real, intent(in) :: x
                 end function func
              end interface

              if (present(loopnum)) then
                 l = loopnum
              else
                 l = 1
              end if

              if (func(a)*func(b)>0.0) then
                 write(*,*)"math::m_bisection: root not in interval"
                 write(*,*)"a:",a,", b:",b,", func(a):",func(a),", func(b):",func(b)
                 write(*,*)"loopnum: ",l
                 stop
              else if (l > maxloop .or. abs(func(a)-func(b))<epsilon) then
                 root = (a+b)/2
              else if (func((a+b)/2.)*func(b)>0) then
                 root = m_bisection(func,a,(a+b)/2.,l+1)
              else
                 root = m_bisection(func,(a+b)/2.,b,l+1)
              end if
       end function m_bisection

       real function m_integral(f, a, b)
              real, intent(in) :: a, b
              real, dimension(:), allocatable :: x, y
              integer :: i

              interface
                 real function f(x)
                        implicit none
                        real, intent(in) :: x
                 end function f
              end interface

              allocate(x(maxloop), y(maxloop))

              do i=1,maxloop
                 x(i) = a + (i-1.)/(maxloop-1.)*(b-a)
                 y(i) = f(x(i))
              end do

              m_integral = m_trapz(x, y)

              deallocate(x, y)
       end function m_integral

       real function m_trapz(x, y, lower, upper)
              real, dimension(:) :: x,y
              integer, optional :: lower, upper
              integer :: l, u, i

              if (present(lower)) then
                 l = lower
              else
                 l = 1
              end if

              if (present(upper)) then
                 u = upper
              else
                 u = size(x)
              end if

              m_trapz = 0.

              do i = l,u-1
                 m_trapz = m_trapz + (x(i+1)-x(i))*(y(i+1)+y(i))/2.0
              end do
       end function m_trapz

       real function m_trapz_poisson_error(x, y, lower, upper)	
              real, dimension(:) :: x,y
              integer, optional :: lower, upper
              integer :: l, u, i

              if (present(lower)) then
                 l = lower
              else
                 l = 1
              end if

              if (present(upper)) then
                 u = upper
              else
                 u = size(x)
              end if

              m_trapz_poisson_error = 0

              do i = l, u-1
                 m_trapz_poisson_error = m_trapz_poisson_error &
                      + (x(i+1)-x(i))/2.0*(sqrt(abs(y(i+1)))+sqrt(abs(y(i))))
              end do
       end function m_trapz_poisson_error

       real function m_sd(v)
              real, dimension(:) :: v
              real :: avg_sum_squared
              integer :: i

              avg_sum_squared = sum(v)**2/size(v)
              ! print*,avg_sum_squared
              m_sd = 0.
              do i=1,size(v)
                 m_sd = m_sd + (v(i)**2)
                 ! print*, m_sd
              end do
              m_sd = sqrt(abs((m_sd-avg_sum_squared)/(size(v)-1.)))
              ! print*, m_sd
       end function m_sd
end module math