* PACKAGE UBISEC  !" $BJ,3dK!$K$h$kJ}Dx<0$N2r$N?dDj(B
*
*"[HIS]  93/07/03(takepiro) Numerical Recipies Sec.9.3 $B$h$j0\?"(B(ZBRENT) 
*
*"       3 $BE@$N4X?tCM$+$iJ}Dx<0$N2r$r?dDj$9$k(B. $B4X?tCM$r7W;;$7K\%k!<%A%s$r(B
*"       $B7+$jJV$7E,MQ$9$k$3$H$K$h$jM?$($i$l$?@:EY$N2r$r7hDj$9$k(B. 
*
**********************************************************************
      SUBROUTINE UBISEC                  !" 2 $BJ,K!$K$h$kJ}Dx<02r$N?dDj(B
     M              ( X1   , X2   , X   , 
     M                FX1  , FX2  , FX  , 
     O                OCNV ,
     I                TOL                  )
*
*" [MODIFY]
      REAL    X1          !" $B2r$NHO0O$=$N(B 1
      REAL    X2          !" $B2r$NHO0O$=$N(B 2
      REAL    X           !" $B?dDj2r(B
*
      REAL    FX1         !" X1 $B$K$*$1$k4X?t$NCM(B
      REAL    FX2         !" X2 $B$K$*$1$k4X?t$NCM(B
      REAL    FX          !" X  $B$K$*$1$k4X?t$NCM(B
*
*" [OUTPUT]
      LOGICAL OCNV        !" $B<}B+$7$?$+H]$+(B
*
*" [INPUT]
      REAL   TOL          !" $B<}B+H=Dj>r7o(B : X $B$N@dBP8m:9(B
*
*" [INTERNAL WORK]
      REAL    XM, TOL1
      REAL    EPS          !" $BIbF0>.?t@:EY(B
      LOGICAL OFIRST
      DATA    OFIRST / .TRUE. /
*
      REAL    R , RP
      REAL*8  DELTPQ
      REAL*8  DELTPD
      REAL    DELTPS
      DATA    DELTPQ / 1.0 D-28 /         !' $B#4G\@:EY$N<}B+>r7o(B
      DATA    DELTPD / 1.0 D-14 /         !'   $BG\@:EY$N<}B+>r7o(B
      DATA    DELTPS / 1.0 E-7  /         !'   $BC1@:EY$N<}B+>r7o(B
*
      SAVE
*
*"  < 0. $BIbF0>.?tE@@:EY%A%'%C%/(B >
*
      IF( OFIRST )THEN
         OFIRST = .FALSE.
         WRITE ( 6, * ) ' BISECTION METHOD : 93/07/03'
*
         R  = 1.
         RP = R + DELTPQ
         IF ( RP .GT. R ) THEN
            WRITE ( 6,* ) ' precision check .. ', RP
            WRITE ( 6,* ) ' REAL*16 is used in UBISEC'
            EPS = DELTPQ
         ELSE
            R  = 1.
            RP = R + DELTPD
            WRITE ( 6,* ) ' precision check .. ', RP
            IF ( RP .GT. R ) THEN
               WRITE ( 6,* ) ' REAL*8 is used in UBISEC'
               EPS = DELTPD
            ELSE
               WRITE ( 6,* ) ' REAL*4 is used in UBISEC'
               EPS = DELTPS
            ENDIF
         ENDIF
      ENDIF
*
*"  < 1. $B<}B+H=Dj(B >
*
      TOL1 = 2.  * EPS * ABS(X) + 0.5 * TOL
      XM   = 0.5 * ( X2 - X1 )
*
      IF(  ABS(XM) .LE. TOL1 .OR. FX. EQ.0. )THEN
         OCNV = .TRUE.
         RETURN
      ENDIF
*
*"  < 2. $BG[CV$,$((B >      
*
      IF(  FX * FX2 .GT. 0. )THEN
         X2  = X
         FX2 = FX
      ELSE
         X1  = X
         FX1 = FX
      ENDIF
*
*"  < 2. $B?7$7$$2r$N?dDj(B >
*
      X = ( X1 + X2 ) / 2.
      OCNV = .FALSE.
*
      RETURN
      END
**********************************************************************
      SUBROUTINE URFLSP            !" $B$O$5$_$&$AK!$K$h$kJ}Dx<02r$N?dDj(B
     M              ( X1   , X2   , X   , 
     M                FX1  , FX2  , FX  , 
     O                OCNV ,
     I                TOL                  )
*
*" [MODIFY]
      REAL    X1          !" $B2r$NHO0O$=$N(B 1
      REAL    X2          !" $B2r$NHO0O$=$N(B 2
      REAL    X           !" $B?dDj2r(B
*
      REAL    FX1         !" X1 $B$K$*$1$k4X?t$NCM(B
      REAL    FX2         !" X2 $B$K$*$1$k4X?t$NCM(B
      REAL    FX          !" X  $B$K$*$1$k4X?t$NCM(B
*
*" [OUTPUT]
      LOGICAL OCNV        !" $B<}B+$7$?$+H]$+(B
*
*" [INPUT]
      REAL    TOL         !" $B<}B+H=Dj>r7o(B : X $B$N@dBP8m:9(B
*
*" [INTERNAL WORK]
      REAL    XM, TOL1
      REAL    EPS         !" $BIbF0>.?t@:EY(B
      LOGICAL OFIRST
      DATA    OFIRST / .TRUE. /
*
      REAL    R , RP
      REAL*8  DELTPQ
      REAL*8  DELTPD
      REAL    DELTPS
      DATA    DELTPQ / 1.0 D-28 /         !' $B#4G\@:EY$N<}B+>r7o(B
      DATA    DELTPD / 1.0 D-14 /         !'   $BG\@:EY$N<}B+>r7o(B
      DATA    DELTPS / 1.0 E-7  /         !'   $BC1@:EY$N<}B+>r7o(B
*
      SAVE
*
*"  < 0. $BIbF0>.?tE@@:EY%A%'%C%/(B >
*
      IF( OFIRST )THEN
         OFIRST = .FALSE.
         WRITE ( 6, * ) ' FALSE POSITION METHOD : 93/07/03'
*
         R  = 1.
         RP = R + DELTPQ
         IF ( RP .GT. R ) THEN
            WRITE ( 6,* ) ' precision check .. ', RP
            WRITE ( 6,* ) ' REAL*16 is used in URFLSP'
            EPS = DELTPQ
         ELSE
            R  = 1.
            RP = R + DELTPD
            WRITE ( 6,* ) ' precision check .. ', RP
            IF ( RP .GT. R ) THEN
               WRITE ( 6,* ) ' REAL*8 is used in URFLSP'
               EPS = DELTPD
            ELSE
               WRITE ( 6,* ) ' REAL*4 is used in URFLSP'
               EPS = DELTPS
            ENDIF
         ENDIF
*
      ENDIF
*
*"  < 1. $B<}B+H=Dj(B >
*
      TOL1 = 2.  * EPS * ABS(X) + 0.5 * TOL
*
      IF ( (X.NE.X1) .AND. (X.NE.X2) )THEN      !" $BE,Ev$J<}B+H=Dj>r7o$+(B?
         IF ( FX * FX1 .GT. 0 ) THEN
            XM   = ABS( X - X1 )
         ELSE
            XM   = ABS( X - X2 )
         ENDIF
*
         IF(  ABS(XM) .LE. TOL1 .OR. FX. EQ.0. )THEN
            OCNV = .TRUE.
            RETURN
         ENDIF
      ENDIF
*
*"  < 2. $BG[CV$,$((B >      
*
      IF(  FX * FX2 .GT. 0. )THEN
         X2  = X
         FX2 = FX
      ELSE
         X1  = X
         FX1 = FX
      ENDIF
*
*"  < 2. $B?7$7$$2r$N?dDj(B >
*
      X = ( FX2 * X1 - FX1 * X2 ) / ( FX2 - FX1 )
      OCNV = .FALSE.
*
      RETURN
      END
**********************************************************************
      SUBROUTINE UBRENT              !" Brent $BK!$K$h$kJ}Dx<0$N2r$N?dDj(B
     M              ( X1   , X2   , X   , 
     M                FX1  , FX2  , FX  , 
     O                OCNV ,
     I                TOL                  )
*
*" [MODIFY]
      REAL   X1           !" $B2r$NHO0O$=$N(B 1
      REAL   X2           !" $B2r$NHO0O$=$N(B 2
      REAL   X            !" $B?dDj2r(B
*
      REAL   FX1          !" X1 $B$K$*$1$k4X?t$NCM(B
      REAL   FX2          !" X2 $B$K$*$1$k4X?t$NCM(B
      REAL   FX           !" X  $B$K$*$1$k4X?t$NCM(B
*
*" [OUTPUT]
      LOGICAL OCNV        !" $B<}B+$7$?$+H]$+(B
*
*" [INPUT]
      REAL   TOL          !" $B<}B+H=Dj>r7o(B
*
*" [INTERNAL WORK]
      REAL   D, E, P, Q, R, S, XM, TOL1
      REAL   EPS          !" $BIbF0>.?t@:EY(B
      LOGICAL OFIRST
      DATA    OFIRST / .TRUE. /
*
      REAL    RR , RP
      REAL*8  DELTPQ
      REAL*8  DELTPD
      REAL    DELTPS
      DATA    DELTPQ / 1.0 D-28 /         !' $B#4G\@:EY$N<}B+>r7o(B
      DATA    DELTPD / 1.0 D-14 /         !'   $BG\@:EY$N<}B+>r7o(B
      DATA    DELTPS / 1.0 E-7  /         !'   $BC1@:EY$N<}B+>r7o(B
*
      SAVE
*
*"  < 0. $BIbF0>.?tE@@:EY%A%'%C%/(B >
*
      IF( OFIRST )THEN
         OFIRST = .FALSE.
         WRITE ( 6, * ) ' BRENT METHOD : 93/07/05'
*
         RR  = 1.
         RP = RR + DELTPQ
         IF ( RP .GT. RR ) THEN
            WRITE ( 6,* ) ' precision check .. ', RP
            WRITE ( 6,* ) ' REAL*16 is used in UBRENT'
            EPS = DELTPQ
         ELSE
            RR  = 1.
            RP = RR + DELTPD
            WRITE ( 6,* ) ' precision check .. ', RP
            IF ( RP .GT. RR ) THEN
               WRITE ( 6,* ) ' REAL*8 is used in UBRENT'
               EPS = DELTPD
            ELSE
               WRITE ( 6,* ) ' REAL*4 is used in UBRENT'
               EPS = DELTPS
            ENDIF
         ENDIF
*
      ENDIF
*
*"  < 1. $BG[CV$,$((B >
*
      IF ( FX * FX2 .GT. 0. ) THEN
         X2  = X1
         FX2 = FX1
         D   = X - X1
         E   = D
      ENDIF
      IF(ABS(FX2).LT.ABS(FX)) THEN
         X1  = X
         X   = X2
         X2  = X1
         FX1 = FX
         FX  = FX2
         FX2 = FX1
      ENDIF
*
*"  < 2. $B<}B+H=Dj(B >
*
      TOL1 = 2.  * EPS * ABS(X) + 0.5*TOL
      XM   = 0.5 * ( X2 - X )
*
      IF(  ABS(XM) .LE. TOL1 .OR. FX. EQ.0. )THEN
         OCNV = .TRUE.
         RETURN
      ENDIF
*
*"  < 3. $B?7$7$$2r$N?dDj(B >
*
      IF(  ABS(E).GE.TOL1 .AND. ABS(FX1).GT.ABS(FX)  ) THEN
         S=  FX / FX1
         IF ( X1 .EQ. X2 ) THEN                  !" $B$O$5$_$&$AK!(B
            P = 2. * XM * S
            Q = 1.-S
         ELSE                                    !" $B5U(B 2 $B<!4X?tJd4V(B
            Q = FX1 / FX2
            R = FX  / FX2
            P = S *( 2.*XM*Q*(Q-R)-(X-X1)*(R-1.) )
            Q = (Q-1.)*(R-1.)*(S-1.)
         ENDIF
         IF(P.GT.0.) Q=-Q
         P=ABS(P)
         IF( 2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q)) ) THEN
            E = D                           !" $B?7$7$$?dDjCM$r<u$1F~$l$k(B
            D = P / Q
         ELSE
            D = XM                          !" 2 $BJ,K!$K$h$k?dDj(B
            E = D
         ENDIF
      ELSE
         D=XM                               !" 2 $BJ,K!$K$h$k?dDj(B
         E=D
      ENDIF
*
      X1  = X                       !"  $B0JA0$N?dDjCM$r6h4V$NC<E@$H$7$FJ]B8(B
      FX1 = FX
*
      IF(  ABS(D) .GT. TOL1 ) THEN  !" $B?7?dDjCM$N7W;;(B
         X  = X + D
      ELSE
         X  = X + SIGN(TOL1,XM)
      ENDIF
*
      OCNV = .FALSE.
*
      RETURN
      END
