!   TITLE:  zero.mac (SOLVE)
!
!               Added "z" prefix to avoid parameter name conflicts.
!               Also return the final y-value, froot = f(root)
!
!   Find the zero of a function using Ridder's method.  Return the
!   result in global parameter "root" and "froot" where froot = f(root).
!   
!   Ref.  Press, "Numerical Recipes in C", 2nd ed., Cambridge, 1992, p358
!
!
!   Usage:  zero,f,x1,x2,offset,xacc
!
!   Arguments:
!       
!          f =  Macro name defining the function to find the root of.
!               The macro must accept the independant variable as 
!               its first argument.  f must return f(x) in global
!               parameter in "yzero".
!
!         x1 =  Lower bracket.  The root must reside between x1 and x2.
!
!         x2 =  Upper bracket.
!
!     offset =  Find where function equals "offset" instead of zero.
!               This will find x=root such that
!
!                   f(root) = offset
!
!       xacc =  The root returned by this macro, "root", will be
!               refined to approximately this accuracy.  Accuracy defaults
!               to 1.0E-4 if xacc is not specified.
!
!
!   Example:    zero,'ztest',.5,2,0,1e-5
!
!               Macro calls of the form
!
!                   ztest,x
!
!               will be repeatedly performed to converge on the root
!               between .5 and 2 with an x accuracy of approximately 1e-5.
!               How close f(x) equals offset will depend on the slope of
!               the curve.
!

*get,prkey_,active,0,prkey
/nopr

zMAXIT_=30
zUNUSED_=-1.11e30

zf_=arg1
zx1_=arg2
zx2_=arg3
zoffset_=arg4
zxacc_=arg5

*if,zoffset_,eq,0,then
    zoffset_=0
*endif

*if,zxacc_,eq,0,then
    zxacc_=1e-4
*endif

zdone_=0                ! Flag set to 1 when converged on root

/out,zero,out
*msg,,zxacc_
Results of zero.mac (xacc = %g):
*msg

*msg,,zoffset_
Loop #___________x________________f(x)___________f(x)-(%g)

*msg

/out

*set,yzero
parsav,,zero,parm
%zf_%,zx1_
parres,change,zero,parm
zfl_=yzero-zoffset_

/out,zero,out,,append
*vwrite,zx1_,yzero,zfl_
(' x1 ',F18.6,E18.4,E18.3)
/out

*set,yzero
parsav,,zero,parm
%zf_%,zx2_
parres,change,zero,parm
zfh_=yzero-zoffset_

/out,zero,out,,append
*vwrite,zx2_,yzero,zfh_
(' x2 ',F18.4,E18.4,E18.3)
/out

and,'zb1_',zfl_>0,zfh_<0
and,'zb2_',zfl_<0,zfh_>0
or,'zb_',zb1_,zb2_

*if,zb_,ne,0,then
    zxl_=zx1_
    zxh_=zx2_
    zans_=zUNUSED_

    *do,zj_,1,zMAXIT_
        zxm_=.5*(zxl_+zxh_)

        *set,yzero
        parsav,,zero,parm
        %zf_%,zxm_
        parres,change,zero,parm
        zfm_=yzero-zoffset_

/out,zero,out,,append
*vwrite,zj_,zxm_,yzero,zfm_
(F3.0,' ',F18.4,E18.4,E18.3)
/out

        zs_=sqrt(zfm_*zfm_-zfl_*zfh_)

        *if,zs_,eq,0,then
            zdone_=1
            *exit
        *endif

        *if,zfl_,ge,zfh_,then
            zsign_=1.0
        *else
            zsign_=-1.0
        *endif

        zxnew_=zxm_+(zxm_-zxl_)*zsign_*zfm_/zs_

        *if,abs(zxnew_-zans_),le,zxacc_,then
            zdone_=1
            *exit
        *endif

        zans_=zxnew_
        *set,yzero
        parsav,,zero,parm
        %zf_%,zans_
        parres,change,zero,parm
        zfnew_=yzero-zoffset_

/out,zero,out,,append
*vwrite,zxnew_,yzero,zfnew_
('    ',F18.4,E18.4,E18.3)
/out

        *if,zfnew_,eq,0,then
            zdone_=1
            *exit
        *endif

        ! if (SIGN(fm,fnew) != fm) {
        *if,zfnew_,ge,0,then
            zb_=abs(zfm_)
        *else
            zb_=-abs(zfm_)
        *endif

        *if,zb_,ne,zfm_,then
            zxl_=zxm_
            zfl_=zfm_
            zxh_=zans_
            zfh_=zfnew_
        *else
          ! if (SIGN(fl,fnew) != fl) {
            *if,zfnew_,ge,0,then
              zb_=abs(zfl_)
            *else
                zb_=-abs(zfl_)
            *endif

            *if,zb_,ne,zfl_,then
                zxh_=zans_
                zfh_=zfnew_
            *else
              ! if (SIGN(fh,fnew) != fh) {
              *if,zfnew_,ge,0,then
                  zb_=abs(zfh_)
              *else
                  zb_=-abs(zfh_)
              *endif

              *if,zb_,ne,zfh_,then
                  zxl_=zans_
                  zfl_=zfnew_
              *else
                *msg,error
                Error in zero.mac: Not supposed to get here.
              *endif
      *endif
        *endif

        *if,abs(zxh_-zxl_),le,zxacc_,then
            zdone_=1
            *exit
        *endif
    *enddo

    *if,zdone_,eq,1,then
        root=zans_
    *else
        *msg,error,zMAXIT_
        Exceeded maximum iterations in %i in zero.mac
    *endif

*else
    *if,zfl_,eq,0,then
        root=zx1_
        zdone_=1
    *endif

    *if,zfh_,eq,0,then
        root=zx2_
        zdone_=1
    *endif

    *if,zdone_,eq,0,then
        *msg,error
        Root must be bracketed in zero.mac
    *endif
*endif


*if,zdone_,eq,1,then

froot=zfnew_+zoffset_

/out,zero,out,,append
*msg

*msg,,root,froot
 *** f(%g) = %g
*msg

*msg,,zj_
Found root in %i loops
/out
*uilist,zero,out

*endif


*set,zMAXIT_
*set,zUNUSED_
*set,zf_
*set,zx1_
*set,zx2_
*set,zoffset_
*set,zxacc_
*set,zdone_
*set,zfl_
*set,zfh_
*set,zfnew_
*set,zb_
*set,zb1_
*set,zb2_
*set,zans_
*set,zj_
*set,zxm_
*set,zxl_
*set,zxh_
*set,zfm_
*set,zs_
*set,zsign_
*set,zxnew_

*if,prkey_,eq,1,then
    /go
*endif