!$Author: joelevin $
!$Date: 2011-03-24 11:33:26 -0500 (Thu, 24 Mar 2011) $
!$Revision: 11724 $
!$HeadURL: https://svn.weru.ksu.edu/weru/weps1/trunk/weps.src/src/lib_wepp/eqroot.for $

      FUNCTION EQROOT(A,ERR)

!     + + + PURPOSE + + +

!     THIS FUNCTION SOLVES THE FOLLOWING EQUATION FOR U:

!       1 - EXP(-U) = A*U, A POSITIVE, U POSITIVE (UNLESS A=1).

!     NEWTON'S METHOD, WITH SPECIAL APRROXIMATIONS FOR SMALL AND 
!     LARGE U, IS USED.  THE RESULTS APPEAR TO BE ACCURATE TO MACHINE
!     PRECISION (REAL*4) AND REQUIRE AT MOST 2 ITERATIONS.

!     LET F(U) = (1 - EXP(-U))/U.

!     FOR A SMALL, F(1/A) = A - A*EXP(-1/A), AND THE LAST TERM IS 
!     SMALL,  SO 1/A IS AN APPROXIMATE SOLUTION.   FOR A <= .06, 
!     THE RELATIVE ERROR IN F AND IN U APPEARS TO BE LESS THAN 
!     REAL*4 PRECISION.

!     LET W = (3 - SQRT(24A - 15))/2.  THEN A = 1 - W/2 + U**2/6.  
!     FROM THE TAYLOR SERIES EXPANSION: 
!          EXP(-U) = 1 - U + U**2/2 - U**3/6 +R(U) 
!     WE GET F(W)=A - R(W)/W, WHERE 0 < R(W) < W**3/24.  FOR A CLOSE 
!     TO 1, W IS NEAR ZERO AND W IS AN APPROXIMATE SOLUTION.  

!     FOR .999 <= A,  THE RELATIVE ERROR IN F AND U APPEARS TO BE 
!     LESS THEN REAL*4 PRECISION.

!     BETWEEN .06 AND .999, NEWTON'S METHOD IS USED,  USING 
!     STARTING VALUES.

!     LET F(U1)=A AND U NEAR U1.    
!      |(A - F(U))/A|  =  |(F(U1) - F(U))/A|  
!         = |F'(C)*U1/A|*|(U1 - U)/U1|,  FOR SOME C BETWEEN U AND U1.
!     SINCE U IS NEAR U1, C AND U1 ARE APPROXIMATED BY U.  THUS FOR 
!     R=|A/(F'(U)*U)|, WE HAVE |(U - U1)/U1| ~ R*|(A - F(U))/A|.  
!     THE RELATIVE ERROR IN F AND U ARE SMALL IF THE RELATIVE ERROR 
!     IN F AND R TIMES THIS ERROR ARE SMALL.
!     EVALUATING F'(U) GIVES R = A/((U - 1)*F(U) - 1).

!     BECAUSE OF THE APPROXIMATIONS IN THE ABOVE, WE USE 1/2 THE
!     SMALLEST REAL*4 NUMBER IN TESTING FOR CONVERGENCE.

!     CALLED FROM DBLEX
!     AUTHOR(S): D. FLANAGAN, J. ASCOUGH
!     VERSION: THIS MODULE TAKEN FROM ASCOUGH STANDALONE IRS CODE
!     DATE CODED:  3-30-2005
!     CODED BY: D. FLANAGAN

!     + + + ARGUMENT DECLARATIONS + + +

      REAL EQROOT, A
      INTEGER ERR  

!     + + + ARGUMENT DEFINITIONS + + +

!     EQROOT   - R*4  RETURNED POSITIVE SOLUTION TO THE EQUATION.
!     A        - R*4  POSITIVE CONSTANT IN THE EQUATION.
!     ERR      - I*2  0: EQUATION SOLVED.
!                     1: NO SOLUTION FOR GIVEN A.

!     + + + LOCAL VARIABLES + + +

      REAL*8 AA, D, E, F, R, S, U

!     + + + END SPECIFICATIONS + + +

!     CHECK A TO SEE IF THERE IS A SOLUTION.

      AA = DBLE(A)
      IF ((AA.LE.0.D0).OR.(1.D0.LT.AA)) THEN
         ERR = 1
         RETURN
      END IF

!     SPECIAL CASE: A=1 (EXACT LIMITING SOLUTION).

      IF (AA.EQ.1.D0) THEN
         ERR = 0
         EQROOT = 0.
         RETURN
      END IF

!     SPECIAL CASE: A SMALL (ANSWER GOOD TO MACHINE PERCISION).

      IF (AA.LE..06D0) THEN
         ERR = 0
         U = 1.D0 / AA
         EQROOT = SNGL(U)
         RETURN
      END IF

!     SPECIAL CASE: A CLOSE TO 1 (ANSWER GOOD TO ABOUT 10 PLACES).

      IF (.999D0.LE.AA) THEN
         ERR = 0
         U = (3.D0/2.D0) - SQRT(6.D0*AA-(15.D0/4.D0))
         EQROOT = SNGL(U)
         RETURN
      END IF

!     ESTIMATE STARTING VALUE FOR U.

      IF (AA.LE..2D0) THEN
         U = 1.D0 / AA
      ELSE IF (AA.LE..5D0) THEN
         U = .968732D0 / AA - 1.55098D0 * AA + .431653D0
      ELSE IF (AA.LE..94D0) THEN
         U = 1.13243D0 / AA - .928240D0 * AA - .207111
      ELSE
         U = (3.D0/2.D0) - SQRT(6.D0*AA-(15.D0/4.D0))
      END IF

!     ITERATE.

   10 CONTINUE
      E = EXP(-U)
      F = (1.D0-E) / U
      D = AA - F
      R = AA / ((U+1.D0)*F-1.D0)
      IF (R.LE.1.D0) THEN
         S = ABS(D/AA)
      ELSE
         S = ABS(R*D/AA)
      END IF
      IF (S.GE..59E-6) THEN
         U = U * (1.D0+D/(E-F))
         GO TO 10
      END IF

!     EXIT WITH SOLUTION.

      ERR = 0
      EQROOT = SNGL(U)
      RETURN
      END