!*****************************************************************************! !*****************************************************************************! !** Complex Root Finding Routine **! !*****************************************************************************! !*****************************************************************************! ! !------------------------------------------------------------------------------ ! Copyright, 2005 ! Greg Howes !------------------------------------------------------------------------------ ! !------------------------------------------------------------------------------ !> FIXME : Add documentation module gk_complex_root implicit none private public :: rtsec contains !------------------------------------------------------------------------------ ! Greg Howes, 2005 !------------------------------------------------------------------------------ ! NOTE: This routine was adapted from f77 routine by Eliot Quataert !> FIXME : Add documentation complex function rtsec(func,x1,x2,xacc,iflag) implicit none integer, parameter :: maxit=1000 interface complex function func(c_) complex, intent(in) :: c_ end function func end interface complex, intent(in) :: x1, x2 real, intent(in) :: xacc integer, intent(out), optional :: iflag complex :: fl, f, swap, dx, maxr, minr complex :: xl integer :: j logical :: limits=.false. !Set limits for solution if (limits) then maxr=x2*10. minr=x1*0.1 endif fl=func(x1) f=func(x2) if(abs(fl).lt.abs(f))then rtsec=x1 xl=x2 swap=fl fl=f f=swap else xl=x1 rtsec=x2 endif do j=1,maxit iflag = j ! write(*,'(a,i4)')'Iteration ',j if (abs(f-fl) .gt. 1.0E-37) then dx=(xl-rtsec)*f/(f-fl) else dx = (x2-x1)/25.0 end if xl=rtsec fl=f ! if (Real(rtsec + dx) .gt. 0.075 .or. Real(rtsec+dx) .lt. 0) then ! rtsec = (x1 + x2)/2.0 ! else !NOTE: Reduce jump by 0.5 to improve convergence GH ! rtsec=rtsec+dx rtsec=rtsec+dx/2. ! end if !Enforce limits if (limits) then !Real limit if (real(rtsec) .gt. real(maxr) .or. real(rtsec) .lt. real(minr)) & rtsec = cmplx(sqrt(real(minr)*real(maxr)), Aimag(rtsec)) !NOTE: aimag(rtsec) should be negative! if (aimag(rtsec) .lt. aimag(maxr) .or. aimag(rtsec) .gt. aimag(minr)) & rtsec = cmplx(real(rtsec), -sqrt(aimag(minr)*aimag(maxr))) endif !LIMIT: Im(rtsec) < 0. ! if (Aimag(rtsec) .gt. 0) then ! rtsec = cmplx(Real(rtsec), -Aimag(rtsec)) ! end if !LIMIT: Re(rtsec) > 0. ! if (Real(rtsec) .lt. 0) then ! rtsec = cmplx(-Real(rtsec), Aimag(rtsec)) ! end if f=func(rtsec) if(abs(dx).lt.xacc.or. Abs(f) .eq. 0.) then ! write(*,'(a,i4)')'Subroutine rtsec Converged on Iteration ',j return endif enddo ! stop ! write(*,'(a,i4)')'Subroutine rtsec stopped at iteration ',j return ! pause 'rtsec exceed maximum iterations' end function rtsec !------------------------------------------------------------------------------ end module gk_complex_root