!+
! Module func_destroyer_mod
!
! Module to get around the NR implementation strategy where "func" (or some other named routine) is used in the code
! where "func" is to be supplied by the user (that is, external to NR), and "func" is *not* passed in as an argument.
!
! This causes problems with putting NR into a shared library since "func" will be an undefined reference.
! The solution is to have modified versions of mnbrak, brent_extended, and dbrent_extended with an extra "func" argument.
!-

module func_destroyer_mod

contains

!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------

	SUBROUTINE mnbrak_extended(ax,bx,cx,fa,fb,fc,func, func_arg)
	USE nrtype; USE nrutil, ONLY : swap
	IMPLICIT NONE
	REAL(dp), INTENT(INOUT) :: ax,bx
	REAL(dp), INTENT(OUT) :: cx,fa,fb,fc

	INTERFACE
		FUNCTION func(x, func_arg)
		USE nrtype
		IMPLICIT NONE
		REAL(dp), INTENT(IN) :: x
	    INTERFACE
		    FUNCTION func_arg(x)
	    	USE nrtype
		    REAL(dp), DIMENSION(:), INTENT(IN) :: x
		    REAL(dp) :: func_arg
		    END FUNCTION func_arg
	    END INTERFACE
		REAL(dp) :: func
		END FUNCTION func

		FUNCTION func_arg(x)
		USE nrtype
		REAL(dp), DIMENSION(:), INTENT(IN) :: x
		REAL(dp) :: func_arg
		END FUNCTION func_arg
	END INTERFACE

	REAL(dp), PARAMETER :: GOLD=1.618034_dp,GLIMIT=100.0_dp,TINY=1.0e-20_dp
	REAL(dp) :: fu,q,r,u,ulim
	fa=func(ax, func_arg)
	fb=func(bx, func_arg)
	if (fb > fa) then
		call swap(ax,bx)
		call swap(fa,fb)
	end if
	cx=bx+GOLD*(bx-ax)
	fc=func(cx, func_arg)
	do
		if (fb < fc) RETURN
		r=(bx-ax)*(fb-fc)
		q=(bx-cx)*(fb-fa)
		u=bx-((bx-cx)*q-(bx-ax)*r)/(2.0_dp*sign(max(abs(q-r),TINY),q-r))
		ulim=bx+GLIMIT*(cx-bx)
		if ((bx-u)*(u-cx) > 0.0) then
			fu=func(u, func_arg)
			if (fu < fc) then
				ax=bx
				fa=fb
				bx=u
				fb=fu
				RETURN
			else if (fu > fb) then
				cx=u
				fc=fu
				RETURN
			end if
			u=cx+GOLD*(cx-bx)
			fu=func(u, func_arg)
		else if ((cx-u)*(u-ulim) > 0.0) then
			fu=func(u, func_arg)
			if (fu < fc) then
				bx=cx
				cx=u
				u=cx+GOLD*(cx-bx)
				call shft(fb,fc,fu,func(u, func_arg))
			end if
		else if ((u-ulim)*(ulim-cx) >= 0.0) then
			u=ulim
			fu=func(u, func_arg)
		else
			u=cx+GOLD*(cx-bx)
			fu=func(u, func_arg)
		end if
		call shft(ax,bx,cx,u)
		call shft(fa,fb,fc,fu)
	end do
	CONTAINS
!BL
	SUBROUTINE shft(a,b,c,d)
	REAL(dp), INTENT(OUT) :: a
	REAL(dp), INTENT(INOUT) :: b,c
	REAL(dp), INTENT(IN) :: d
	a=b
	b=c
	c=d
	END SUBROUTINE shft
	END SUBROUTINE mnbrak_extended

!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------

	FUNCTION brent_extended(ax,bx,cx,func, func_arg, tol,xmin)
	USE nrtype; USE nrutil, ONLY : nrerror
	IMPLICIT NONE
	REAL(dp), INTENT(IN) :: ax,bx,cx,tol
	REAL(dp), INTENT(OUT) :: xmin
	REAL(dp) :: brent_extended

	INTERFACE
		FUNCTION func(x, func_arg)
		USE nrtype
		IMPLICIT NONE
		REAL(dp), INTENT(IN) :: x
	    INTERFACE
		    FUNCTION func_arg(x)
	    	USE nrtype
		    REAL(dp), DIMENSION(:), INTENT(IN) :: x
		    REAL(dp) :: func_arg
		    END FUNCTION func_arg
	    END INTERFACE
		REAL(dp) :: func
		END FUNCTION func

		FUNCTION func_arg(x)
		USE nrtype
		REAL(dp), DIMENSION(:), INTENT(IN) :: x
		REAL(dp) :: func_arg
		END FUNCTION func_arg
	END INTERFACE

	INTEGER(I4B), PARAMETER :: ITMAX=100
	REAL(dp), PARAMETER :: CGOLD=0.3819660_dp,ZEPS=1.0e-3_dp*epsilon(ax)
	INTEGER(I4B) :: iter
	REAL(dp) :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
	a=min(ax,cx)
	b=max(ax,cx)
	v=bx
	w=v
	x=v
	e=0.0
	fx=func(x, func_arg)
	fv=fx
	fw=fx
	do iter=1,ITMAX
		xm=0.5_dp*(a+b)
		tol1=tol*abs(x)+ZEPS
		tol2=2.0_dp*tol1
		if (abs(x-xm) <= (tol2-0.5_dp*(b-a))) then
			xmin=x
			brent_extended=fx
			RETURN
		end if
		if (abs(e) > tol1) then
			r=(x-w)*(fx-fv)
			q=(x-v)*(fx-fw)
			p=(x-v)*q-(x-w)*r
			q=2.0_dp*(q-r)
			if (q > 0.0) p=-p
			q=abs(q)
			etemp=e
			e=d
			if (abs(p) >= abs(0.5_dp*q*etemp) .or. &
				p <= q*(a-x) .or. p >= q*(b-x)) then
				e=merge(a-x,b-x, x >= xm )
				d=CGOLD*e
			else
				d=p/q
				u=x+d
				if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
			end if
		else
			e=merge(a-x,b-x, x >= xm )
			d=CGOLD*e
		end if
		u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
		fu=func(u, func_arg)
		if (fu <= fx) then
			if (u >= x) then
				a=x
			else
				b=x
			end if
			call shft(v,w,x,u)
			call shft(fv,fw,fx,fu)
		else
			if (u < x) then
				a=u
			else
				b=u
			end if
			if (fu <= fw .or. w == x) then
				v=w
				fv=fw
				w=u
				fw=fu
			else if (fu <= fv .or. v == x .or. v == w) then
				v=u
				fv=fu
			end if
		end if
	end do
	call nrerror('brent_extended: exceed maximum iterations')
	CONTAINS
!BL
	SUBROUTINE shft(a,b,c,d)
	REAL(dp), INTENT(OUT) :: a
	REAL(dp), INTENT(INOUT) :: b,c
	REAL(dp), INTENT(IN) :: d
	a=b
	b=c
	c=d
	END SUBROUTINE shft
	END FUNCTION brent_extended

!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------

	FUNCTION dbrent_extended(ax,bx,cx,func, func_arg, dfunc, dfunc_arg, tol,xmin)
	USE nrtype; USE nrutil, ONLY : nrerror
	IMPLICIT NONE
	REAL(dp), INTENT(IN) :: ax,bx,cx,tol
	REAL(dp), INTENT(OUT) :: xmin
	REAL(dp) :: dbrent_extended

	INTERFACE
		FUNCTION func(x, func_arg)
		USE nrtype
		IMPLICIT NONE
		REAL(dp), INTENT(IN) :: x
	    INTERFACE
		    FUNCTION func_arg(x)
	    	USE nrtype
		    REAL(dp), DIMENSION(:), INTENT(IN) :: x
		    REAL(dp) :: func_arg
		    END FUNCTION func_arg
	    END INTERFACE
		REAL(dp) :: func
		END FUNCTION func
!BL
   		FUNCTION func_arg(x)
   		USE nrtype
   		REAL(dp), DIMENSION(:), INTENT(IN) :: x
   		REAL(dp) :: func_arg
   		END FUNCTION func_arg

		FUNCTION dfunc(x, dfunc_arg)
		USE nrtype
		IMPLICIT NONE
		REAL(dp), INTENT(IN) :: x
    	INTERFACE
    		FUNCTION dfunc_arg(x)
    		USE nrtype
    		REAL(dp), DIMENSION(:), INTENT(IN) :: x
    		REAL(dp), DIMENSION(size(x)) :: dfunc_arg
    		END FUNCTION dfunc_arg
    	END INTERFACE
		REAL(dp) :: dfunc
		END FUNCTION dfunc

		FUNCTION dfunc_arg(x)
		USE nrtype
		REAL(dp), DIMENSION(:), INTENT(IN) :: x
		REAL(dp), DIMENSION(size(x)) :: dfunc_arg
		END FUNCTION dfunc_arg
	END INTERFACE

	INTEGER(I4B), PARAMETER :: ITMAX=100
	REAL(dp), PARAMETER :: ZEPS=1.0e-3_dp*epsilon(ax)
	INTEGER(I4B) :: iter
	REAL(dp) :: a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,&
		u,u1,u2,v,w,x,xm
	LOGICAL :: ok1,ok2
	a=min(ax,cx)
	b=max(ax,cx)
	v=bx
	w=v
	x=v
	e=0.0
	fx=func(x, func_arg)
	fv=fx
	fw=fx
	dx=dfunc(x, dfunc_arg)
	dv=dx
	dw=dx
	do iter=1,ITMAX
		xm=0.5_dp*(a+b)
		tol1=tol*abs(x)+ZEPS
		tol2=2.0_dp*tol1
		if (abs(x-xm) <= (tol2-0.5_dp*(b-a))) exit
		if (abs(e) > tol1) then
			d1=2.0_dp*(b-a)
			d2=d1
			if (dw /= dx) d1=(w-x)*dx/(dx-dw)
			if (dv /= dx) d2=(v-x)*dx/(dx-dv)
			u1=x+d1
			u2=x+d2
			ok1=((a-u1)*(u1-b) > 0.0) .and. (dx*d1 <= 0.0)
			ok2=((a-u2)*(u2-b) > 0.0) .and. (dx*d2 <= 0.0)
			olde=e
			e=d
			if (ok1 .or. ok2) then
				if (ok1 .and. ok2) then
					d=merge(d1,d2, abs(d1) < abs(d2))
				else
					d=merge(d1,d2,ok1)
				end if
				if (abs(d) <= abs(0.5_dp*olde)) then
					u=x+d
					if (u-a < tol2 .or. b-u < tol2) &
						d=sign(tol1,xm-x)
				else
					e=merge(a,b, dx >= 0.0)-x
					d=0.5_dp*e
				end if
			else
				e=merge(a,b, dx >= 0.0)-x
				d=0.5_dp*e
			end if
		else
			e=merge(a,b, dx >= 0.0)-x
			d=0.5_dp*e
		end if
		if (abs(d) >= tol1) then
			u=x+d
			fu=func(u, func_arg)
		else
			u=x+sign(tol1,d)
			fu=func(u, func_arg)
			if (fu > fx) exit
		end if
		du=dfunc(u, dfunc_arg)
		if (fu <= fx) then
			if (u >= x) then
				a=x
			else
				b=x
			end if
			call mov3(v,fv,dv,w,fw,dw)
			call mov3(w,fw,dw,x,fx,dx)
			call mov3(x,fx,dx,u,fu,du)
		else
			if (u < x) then
				a=u
			else
				b=u
			end if
			if (fu <= fw .or. w == x) then
				call mov3(v,fv,dv,w,fw,dw)
				call mov3(w,fw,dw,u,fu,du)
			else if (fu <= fv .or. v == x .or. v == w) then
				call mov3(v,fv,dv,u,fu,du)
			end if
		end if
	end do
	if (iter > ITMAX) call nrerror('dbrent_extended: exceeded maximum iterations')
	xmin=x
	dbrent_extended=fx
	CONTAINS
!BL
	SUBROUTINE mov3(a,b,c,d,e,f)
	REAL(dp), INTENT(IN) :: d,e,f
	REAL(dp), INTENT(OUT) :: a,b,c
	a=d
	b=e
	c=f
	END SUBROUTINE mov3
	END FUNCTION dbrent_extended

!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------

	SUBROUTINE lnsrch_extended(xold,fold,g,p,x,f,stpmax,check,func, func_arg)
	USE nrtype; USE nrutil, ONLY : assert_eq,nrerror,vabs
	IMPLICIT NONE
	REAL(dp), DIMENSION(:), INTENT(IN) :: xold,g
	REAL(dp), DIMENSION(:), INTENT(INOUT) :: p
	REAL(dp), INTENT(IN) :: fold,stpmax
	REAL(dp), DIMENSION(:), INTENT(OUT) :: x
	REAL(dp), INTENT(OUT) :: f
	LOGICAL(LGT), INTENT(OUT) :: check
	INTERFACE
		FUNCTION func(x, func_arg)
		USE nrtype
		IMPLICIT NONE
	    INTERFACE
	    	FUNCTION func_arg(x)
	    	USE nrtype
    		IMPLICIT NONE
	    	REAL(dp), DIMENSION(:), INTENT(IN) :: x
	    	REAL(dp), DIMENSION(size(x)) :: func_arg
	    	END FUNCTION func_arg
    	END INTERFACE
		REAL(dp) :: func
		REAL(dp), DIMENSION(:), INTENT(IN) :: x
		END FUNCTION func

		FUNCTION func_arg(x)
		USE nrtype
		IMPLICIT NONE
		REAL(dp), DIMENSION(:), INTENT(IN) :: x
		REAL(dp), DIMENSION(size(x)) :: func_arg
		END FUNCTION func_arg
	END INTERFACE

	REAL(dp), PARAMETER :: ALF=1.0e-4_dp,TOLX=epsilon(x)
	INTEGER(I4B) :: ndum
	REAL(dp) :: a,alam,alam2,alamin,b,disc,f2,pabs,rhs1,rhs2,slope,tmplam
	ndum=assert_eq(size(g),size(p),size(x),size(xold),'lnsrch')
	check=.false.
	pabs=vabs(p(:))
	if (pabs > stpmax) p(:)=p(:)*stpmax/pabs
	slope=dot_product(g,p)
	if (slope >= 0.0) call nrerror('roundoff problem in lnsrch')
	alamin=TOLX/maxval(abs(p(:))/max(abs(xold(:)),1.0_dp))
	alam=1.0
	do
		x(:)=xold(:)+alam*p(:)
		f=func(x, func_arg)
		if (alam < alamin) then
			x(:)=xold(:)
			check=.true.
			RETURN
		else if (f <= fold+ALF*alam*slope) then
			RETURN
		else
			if (alam == 1.0) then
				tmplam=-slope/(2.0_dp*(f-fold-slope))
			else
				rhs1=f-fold-alam*slope
				rhs2=f2-fold-alam2*slope
				a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2)
				b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/&
					(alam-alam2)
				if (a == 0.0) then
					tmplam=-slope/(2.0_dp*b)
				else
					disc=b*b-3.0_dp*a*slope
					if (disc < 0.0) then
						tmplam=0.5_dp*alam
					else if (b <= 0.0) then
						tmplam=(-b+sqrt(disc))/(3.0_dp*a)
					else
						tmplam=-slope/(b+sqrt(disc))
					end if
				end if
				if (tmplam > 0.5_dp*alam) tmplam=0.5_dp*alam
			end if
		end if
		alam2=alam
		f2=f
		alam=max(tmplam,0.1_dp*alam)
	end do
	END SUBROUTINE lnsrch_extended

end module