Module Root

! procedury pro hledani korene rovnice f(x)=0
! v intervalu <x1,x2> s chybou xacc
! v maxit je na vstupu maximalni pocet povolenych iteraci a
!            na vystupu pocet skutecne provedenych iteraci
! Vystupni promenna ierr nese informaci o prubehu vypoctu
! ierr = 0 , vypocet bez problemu
! ierr = 1 , byl dosazen maximalni pocet iteraci zadany v maxit
! ierr = 2 , a) funkce nema na krajich intervalu ruzna znamenka
!               pro RTBisect
!            b) vysla ze zadaneho intervalu pro RTNewton
! ierr = 3 , jen pro RTNewton: kombinace ierr=1 + ierr=2
!
! Procedury vyzaduji deklaraci funkce F(x),
! RTNewton navic DF(x), ktera vraci dF/dx
! 
! Pomocna procedura Findb muze najit interval, ktery obsahuje koren


use std_type

implicit none
private
public :: Findb, RTBisect, RTSecant, RTNewton


contains

subroutine Findb(h,a,b,F)
	! Najde interval <a,b> v nemz ma F koren, takze F(a)*F(b)<0
	! Vstup : h = pocatecni krok pro hledani
	!         a = starovaci bod od ktereho zacina hledat
	! Vystup: hledany interval <a,b>
	!         h = posledni hodnota h	
	real(kind=DP),intent(inout) :: a,h
	real(kind=DP),intent(out) :: b
	interface
		function F(x) result(v)
		use std_type
		real(kind=DP),intent(in) :: x
		real(kind=DP) :: v
		end function F
	end interface
	
	real(kind=DP) :: fa,fb
	fa = F(a)
	b = a+h
	fb = F(b)
	do
		if (fa*fb<0) exit
		b = b+h
		fb = F(b)
		if ((fa*fb)>0) then
			if (abs(fa)<abs(fb)) then
				h = -2*h
			else
				h = 2*h
				a = b
				fa =fb
			end if
		end if
	end do
end subroutine Findb


subroutine RTBisect(x1,x2,xacc,koren,ierr,F,maxit)
	! Metoda bisekce
	real(kind=DP),intent(in) :: x1,x2,xacc
	real(kind=DP),intent(out) :: koren
	integer,intent(out) :: ierr
	integer,intent(inout) :: maxit
	interface
		function F(x) result(v)
		use std_type
		real(kind=DP),intent(in) :: x
		real(kind=DP) :: v
		end function F
	end interface
	
	real(kind=DP) :: dx,fh,fmid,xmid,rtb
	integer :: j
	
	ierr = 0
	fmid = F(x2)
	fh = F(x1)
	if (fh*fmid>=0) then
		ierr = 2
		koren = 0
	else   ! interval uzavira koren
		if (fh<0) then
			rtb = x1
			dx = x2-x1
		else
			rtb = x2
			dx = x1-x2
		end if
		j = 0
		do
			j = j+1
			dx = dx/2
			xmid = rtb+dx
			fmid = F(xmid)
			if (fmid<0) rtb = xmid
			if ((abs(dx)<xacc) .or. (fmid==0) .or. (j==maxit)) exit
		end do
		if (j==maxit) ierr = 1
		maxit = j
		koren = rtb
	end if
end subroutine RTBisect



subroutine RTSecant(x1,x2,xacc,koren,ierr,F,maxit)
	! Metoda secen
	real(kind=DP),intent(in) :: x1,x2,xacc
	real(kind=DP),intent(out) :: koren
	integer,intent(out) :: ierr
	integer,intent(inout) :: maxit
	interface
		function F(x) result(v)
		use std_type
		real(kind=DP),intent(in) :: x
		real(kind=DP) :: v
		end function F
	end interface		
	
	real(kind=DP) :: dx,fg,fh,swap,xh,rts
	integer :: j
	
	ierr = 0
	fh = F(x1)
	fg = F(x2)
	if (abs(fh)<abs(fg)) then
		rts = x1
		xh = x2
		swap = fh
		fh = fg
		fg = swap
	else
		xh = x1
		rts = x2
	end if
	j = 0
	do
		j = j+1
		dx = (xh-rts)*fg/(fg-fh)
		xh = rts
		fh = fg
		rts = rts+dx
		fg = F(rts)
		if ((abs(dx)<xacc) .or. (fg==0) .or. (j==maxit)) exit
	end do
	if (j==maxit) ierr = 1
	maxit = j
	koren = rts
end subroutine RTSecant

		
subroutine RTNewton(x1,x2,xacc,koren,ierr,F,DF,maxit)
	! Newtonova metoda
	real(kind=DP),intent(in) :: x1,x2,xacc
	real(kind=DP),intent(out) :: koren
	integer,intent(out) :: ierr
	integer,intent(inout) :: maxit
	interface
		function F(x) result(v)
		use std_type
		real(kind=DP),intent(in) :: x
		real(kind=DP) :: v
		end function F
	end interface
	interface	
		function DF(x) result(v)
		use std_type
		real(kind=DP),intent(in) :: x
		real(kind=DP) :: v
		end function DF
	end interface
	
	real(kind=DP) :: fd,dx,fh,rtn
	integer :: j
	
	ierr = 0
	rtn =(x1+x2)/2
	j = 0
	do
		j = j+1
		fh = F(rtn)
		fd = DF(rtn)
		dx = fh/fd
		rtn = rtn-dx
		if ((x1-rtn)*(rtn-x2)<0) ierr = 2
		if ((abs(dx)<xacc).or.(j==maxit).or.(ierr==2)) exit
	end do
	if (j==maxit) ierr = ierr+1
	maxit = j
	if (ierr==0) then 
		koren = rtn
	else
		koren = 0
	end if
end subroutine RTNewton



end module Root
