CC+------------------------------------------------------------------+CC
CC+------------------------------------------------------------------+CC
CC+                                                                  +CC
CC+              Lawrence Livermore National Laboratory              +CC
CC+                                                                  +CC
CC+     Livermore Computing Center  Mathematical Software Library    +CC
CC+    Mathematical Software Service Library (MSSL) -- Class Three   +CC
CC+                                                                  +CC
CC+------------------------------------------------------------------+CC
CC+                                                                  +CC
CC+  Each Class Three routine is guaranteed only to meet minimal     +CC
CC+  documentation and programming standards.  Although the quality  +CC
CC+  of a particular routine may be high, this is not assured.  A    +CC
CC+  Class Three routine is recommended only in cases where no       +CC
CC+  suitable routine exists in other libraries (e.g., MSSL or       +CC
CC+  SLATEC).                                                        +CC
CC+                                                                  +CC
CC+  These routines are distributed exclusively for use in support   +CC
CC+  of LLNL programs.  Check with the LLNL Code Release Center or   +CC
CC+  the LC Client Services HotLine, (510)422-4531, before moving    +CC
CC+  this source code to a non-LLNL system.                          +CC
CC+                                                                  +CC
CC+  Support for Class Three routines is not guaranteed, but is      +CC
CC+  available in many cases.                                        +CC
CC+                                                                  +CC
CC+  +------------------------------------------------------------+  +CC
CC+  +                        N O T I C E                         +  +CC
CC+  +------------------------------------------------------------+  +CC
CC+  +  This report was prepared as an account of work sponsored  +  +CC
CC+  +  by the United States government.  Neither the United      +  +CC
CC+  +  States government nor any of their employees, nor any of  +  +CC
CC+  +  their contractors, subcontractors, or their employees,    +  +CC
CC+  +  makes any warranty, expressed or implied, or assumes any  +  +CC
CC+  +  legal liability or responsibility for the accuracy,       +  +CC
CC+  +  completeness or usefulness of any information, apparatus, +  +CC
CC+  +  product or process disclosed, or represents that its use  +  +CC
CC+  +  would not infringe privately-owned rights.                +  +CC
CC+  +------------------------------------------------------------+  +CC
CC+                                                                  +CC
CC+  Please refer questions regarding this software to the LC        +CC
CC+  Client Services Hotline, (510)422-4531, or use the NMG "mail"   +CC
CC+  command.                                                        +CC
CC+                                                                  +CC
CC+------------------------------------------------------------------+CC
CC+------------------------------------------------------------------+CC


      real function rfl(x,y,z,ier) 
c***begin prologue  rfl
c***date written   790801   (yymmdd)
c***revision date  831228   (yymmdd)
c***category no.  c14
c***keywords  duplication theorem,elliptic integral,incomplete,complete,
c             integral of the first kind,taylor series
c***author  carlson, b.c., ames laboratory-doe
c             iowa state university, ames, iowa  50011
c           notis, e.m., ames laboratory-doe
c             iowa state university, ames, iowa  50011
c           pexton, r.l., lawrence livermore national laboratory
c             livermore, california  94550
c***purpose  incomplete or complete elliptic integral of the 1st kind.
c            for x, y, and z nonnegative and at most one of them zero,
c            rfl(x,y,z) = integral from zero to infinity of
c                                -1/2     -1/2     -1/2
c                      (1/2)(t+x)    (t+y)    (t+z)    dt,
c            if one of them is zero, the integral is complete.
c***description
c
c   1.     rfl
c          evaluate an incomplete (or complete) elliptic integral
c          of the first kind
c          standard fortran function routine
c          single precision version
c          the routine calculates an approximation result to
c          rfl(x,y,z) = integral from zero to infinity of
c
c                               -1/2     -1/2     -1/2
c                     (1/2)(t+x)    (t+y)    (t+z)    dt,
c
c          where x, y, and z are nonnegative and at most one of them
c          is zero.  if one of them is zero, the integral is complete.
c          the duplication theorem is iterated until the variables are
c          nearly equal, and the function is then expanded in taylor
c          series to fifth order.
c
c   2.     calling sequence
c          rfl( x, y, z, ier )
c
c          parameters on entry
c          values assigned by the calling routine
c
c          x      - single precision, nonnegative variable
c
c          y      - single precision, nonnegative variable
c
c          z      - single precision, nonnegative variable
c
c
c
c          on return     (values assigned by the rfl routine)
c
c          rfl     - single precision approximation to the integral
c
c          ier    - integer
c
c                   ier = 0 normal and reliable termination of the
c                           routine.  it is assumed that the requested
c                           accuracy has been achieved.
c
c                   ier >  0 abnormal termination of the routine
c
c          x, y, z are unaltered.
c
c
c   3.    error messages
c
c         value of ier assigned by the rfl routine
c
c                  value assigned         error message printed
c                  ier = 1                amin1(x,y,z) .lt. 0.0e0
c                      = 2                amin1(x+y,x+z,y+z) .lt. lolim
c                      = 3                amax1(x,y,z) .gt. uplim
c
c
c
c   4.     control parameters
c
c                  values of lolim,uplim,and errtol are set by the
c                  routine.
c
c          lolim and uplim determine the valid range of x, y and z
c
c          lolim  - lower limit of valid arguments
c
c                   not less than 5 * (machine minimum).
c
c          uplim  - upper limit of valid arguments
c
c                   not greater than (machine maximum) / 5.
c
c
c                     acceptable values for:   lolim      uplim
c                     ibm 360/370 series   :   3.0e-78     1.0e+75
c                     cdc 6000/7000 series :   1.0e-292    1.0e+321
c                     univac 1100 series   :   1.0e-37     1.0e+37
c                     cray                 :   2.3e-2466   1.09e+2465
c                     vax 11 series        :   1.5e-38     3.0e+37
c
c
c
c          errtol determines the accuracy of the answer
c
c                 the value assigned by the routine will result
c                 in solution precision within 1-2 decimals of
c                 "machine precision".
c
c
c
c          errtol - relative error due to truncation is less than
c                   errtol ** 6 / (4 * (1-errtol)  .
c
c
c
c              the accuracy of the computed approximation to the inte-
c              gral can be controlled by choosing the value of errtol.
c              truncation of a taylor series after terms of fifth order
c              introduces an error less than the amount shown in the
c              second column of the following table for each value of
c              errtol in the first column.  in addition to the trunca-
c              tion error there will be round-off error, but in prac-
c              tice the total error from both sources is usually less
c              than the amount given in the table.
c
c
c
c
c
c          sample choices:  errtol   relative truncation
c                                    error less than
c                           1.0e-3    3.0e-19
c                           3.0e-3    2.0e-16
c                           1.0e-2    3.0e-13
c                           3.0e-2    2.0e-10
c                           1.0e-1    3.0e-7
c
c
c                    decreasing errtol by a factor of 10 yields six more
c                    decimal digits of accuracy at the expense of one or
c                    two more iterations of the duplication theorem.
c***long description
c
c   rfl special comments
c
c
c
c          check by addition theorem: rfl(x,x+z,x+w) + rfl(y,y+z,y+w)
c          = rfl(0,z,w), where x,y,z,w are positive and x * y = z * w.
c
c
c          on input:
c
c          x, y, and z are the variables in the integral rfl(x,y,z).
c
c
c          on output:
c
c
c          x, y, and z are unaltered.
c
c
c
c          ********************************************************
c
c          warning: changes in the program may improve speed at the
c          expense of robustness.
c
c
c
c   special functions via rfl
c
c
c                  legendre form of elliptic integral of 1st kind
c                  ----------------------------------------------
c
c
c                                             2         2   2
c                  f(phi,k) = sin(phi) rfl(cos (phi),1-k sin (phi),1)
c
c
c                                  2
c                  k(k) = rfl(0,1-k ,1)
c
c
c
c                  bulirsch form of elliptic integral of 1st kind
c                  ----------------------------------------------
c
c
c                                          2 2    2
c                  el1(x,kc) = x rfl(1,1+kc x ,1+x )
c
c
c
c
c                  lemniscate constant a
c                  ---------------------
c
c
c                       1      4 -1/2
c                  a = int (1-s )    ds = rfl(0,1,2) = rfl(0,2,1)
c                       0
c
c
c    -------------------------------------------------------------------
c          subroutines or functions needed
c              - errpex
c              - r1pext
c              - fortran abs,amax1,amin1,sqrt
c***references  carlson, b.c. and notis, e.m.
c                 algorithms for incomplete elliptic integrals
c                 acm transactions on mathematical software,vol.7,no.3,
c                 sept, 1981, pages 398-403
c               carlson, b.c.
c                 computing elliptic integrals by duplication
c                 numer. math. 33, (1979), 1-16
c               carlson, b.c.
c                 elliptic integrals of the first kind
c                 siam j. math. anal. 8 (1977), 231-242
c***routines called  r1pext,errpex
c***end prologue  rfl
      integer ier,itodo
      real lolim, uplim, epslon, errtol
      real c1, c2, c3, e2, e3, lamda
      real mu, s, x, xn, xndev
      real xnroot, y, yn, yndev, ynroot, z, zn, zndev,
     * znroot
c
      character msg*80
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c     the original routine is designed for cft
c
c
c     for  civic  users replace the above character statement
c
c
c        '  character msg*80  '   with the following :
c
c
c           dimension msg(10)
c
c
c     for  chat  users replace the above character statement
c
c
c        '  character msg*80  '   with the following :
c
c
c           dimension msg(8)
c
c
c      other changes are listed subsequently.
c
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c
c
c
      data itodo/1/
c
c
c
c***first executable statement  rfl
      if(itodo.eq.1)then
c
c
      errtol=(4.0*r1pext(3))**(1.0/6.0)
c
c
      lolim = 5.0e0 * (r1pext(1) )
c
      uplim = (r1pext(2) ) / 5.0e0
c
      c1=1.0e0/24.0e0
c
      c2=3.0e0/44.0e0
c
      c3=1.0e0/14.0e0
c
c
      itodo=0
c
      end if
c
c
c         call error handler if necessary.
c
c
    5 rfl=0.0e0
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c     for civic users insert the following three statements:
c
c
c          do 9191 igogo=1,10
c
c          msg(igogo)='        '
c
c9191 continue
c
c
c     for chat users insert the following three statements:
c
c
c          do 9191 igogo=1,8
c
c          msg(igogo)='          '
c
c9191 continue
c
c
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
      if( amin1(x,y,z).lt.0.0e0) then
      ier=1
      msg= 'rfl - error: amin1(x,y,z).lt.0.0e0 (where x=r1, y=r2, z=r3)'
      call errpex (msg,ier,x,y,z,0.0,0.0)
      return
      end if
      if (amin1(x+y,x+z,y+z).lt.lolim) then
      ier=2
      msg='rfl-error:amin1(x+y,x+z,y+z).lt.lolim(where x=r1,y=r2,z=r3)'
      call errpex (msg,ier,x,y,z,0.0,0.0)
      return
      end if
      if (amax1(x,y,z).gt.uplim) then
      ier=3
      msg='rfl - error: amax1(x,y,z).gt.uplim (where x=r1, y=r2, z=r3)'
      call errpex (msg,ier,x,y,z,0.0,0.0)
      return
      end if
c
c
c
   20 ier = 0
      xn = x
      yn = y
      zn = z
c
c
c
   30 mu = (xn+yn+zn)/3.0e0
      xndev = 2.0e0 - (mu+xn)/mu
      yndev = 2.0e0 - (mu+yn)/mu
      zndev = 2.0e0 - (mu+zn)/mu
      epslon = amax1( abs(xndev), abs(yndev), abs(zndev))
      if (epslon.lt.errtol) go to 40
      xnroot =  sqrt(xn)
      ynroot =  sqrt(yn)
      znroot =  sqrt(zn)
      lamda = xnroot*(ynroot+znroot) + ynroot*znroot
      xn = (xn+lamda)*0.250e0
      yn = (yn+lamda)*0.250e0
      zn = (zn+lamda)*0.250e0
      go to 30
c
c
c
   40 e2 = xndev*yndev - zndev*zndev
      e3 = xndev*yndev*zndev
      s = 1.0e0 + (c1*e2-0.10e0-c2*e3)*e2 + c3*e3
      rfl = s/ sqrt(mu)
c
   50 return
c
c
c
      end
      subroutine errpex(messg,ierr,r1,r2,r3,r4,r5)    
c
c
      character messg*80
c
c
c
c***first executable statement  errpex
c
c
      write (6,101) messg
  101 format (a)
c
      write(6,102) ierr,r1,r2,r3,r4,r5
  102 format (2x,i3,3x,2e16.7,/,3x,3e16.7)
      return
      end
      real function r1pext(i)      
c
c
      character msg*80
c
      integer small(2)
      integer large(2)
      integer right(2)
      integer diver(2)
      integer log10(2)
c
      real rmach(5)
c
      equivalence (rmach(1),small(1))
      equivalence (rmach(2),large(1))
      equivalence (rmach(3),right(1))
      equivalence (rmach(4),diver(1))
      equivalence (rmach(5),log10(1))
c
c     machine constants for the ibm-pc
c
       rmach(1)  =2.0**(-149)
       rmach(2)  = 2.0**(127)
       rmach(3)  =2.0**(-23)
       rmach(4)  =2.0**(-22)
       rmach(5)  = 0.30102999566
c
c
c     machine constants for the cdc 6000/7000 series.
c
c     data rmach(1) / 00014000000000000000b /
c     data rmach(2) / 37767777777777777777b /
c     data rmach(3) / 16404000000000000000b /
c     data rmach(4) / 16414000000000000000b /
c     data rmach(5) / 17164642023241175720b /
c
c     machine constants for the cray 1
c
c      data rmach(1) / 200004000000000000000b /
c      data rmach(2) / 577777777777777777777b /
c      data rmach(3) / 377214000000000000000b /
c      data rmach(4) / 377224000000000000000b /
c      data rmach(5) / 377774642023241175720b /
c
c
c
c
      if (i .lt. 1  .or.  i .gt. 5)  then
         msg=  'r mach -- i out of bounds'
         call errpex (msg , i, 0.0, 0.0, 0.0, 0.0, 0.0)
      end if
c
      r1pext = rmach(i)
      return
c
      end

CC+------------------------------------------------------------------+CC
CC+------------------------------------------------------------------+CC
CC+                                                                  +CC
CC+              Lawrence Livermore National Laboratory              +CC
CC+                                                                  +CC
CC+     Livermore Computing Center  Mathematical Software Library    +CC
CC+    Mathematical Software Service Library (MSSL) -- Class Three   +CC
CC+                                                                  +CC
CC+------------------------------------------------------------------+CC
CC+                                                                  +CC
CC+  Each Class Three routine is guaranteed only to meet minimal     +CC
CC+  documentation and programming standards.  Although the quality  +CC
CC+  of a particular routine may be high, this is not assured.  A    +CC
CC+  Class Three routine is recommended only in cases where no       +CC
CC+  suitable routine exists in other libraries (e.g., MSSL or       +CC
CC+  SLATEC).                                                        +CC
CC+                                                                  +CC
CC+  These routines are distributed exclusively for use in support   +CC
CC+  of LLNL programs.  Check with the LLNL Code Release Center or   +CC
CC+  the LC Client Services HotLine, (510)422-4531, before moving    +CC
CC+  this source code to a non-LLNL system.                          +CC
CC+                                                                  +CC
CC+  Support for Class Three routines is not guaranteed, but is      +CC
CC+  available in many cases.                                        +CC
CC+                                                                  +CC
CC+  +------------------------------------------------------------+  +CC
CC+  +                        N O T I C E                         +  +CC
CC+  +------------------------------------------------------------+  +CC
CC+  +  This report was prepared as an account of work sponsored  +  +CC
CC+  +  by the United States government.  Neither the United      +  +CC
CC+  +  States government nor any of their employees, nor any of  +  +CC
CC+  +  their contractors, subcontractors, or their employees,    +  +CC
CC+  +  makes any warranty, expressed or implied, or assumes any  +  +CC
CC+  +  legal liability or responsibility for the accuracy,       +  +CC
CC+  +  completeness or usefulness of any information, apparatus, +  +CC
CC+  +  product or process disclosed, or represents that its use  +  +CC
CC+  +  would not infringe privately-owned rights.                +  +CC
CC+  +------------------------------------------------------------+  +CC
CC+                                                                  +CC
CC+  Please refer questions regarding this software to the LC        +CC
CC+  Client Services Hotline, (510)422-4531, or use the NMG "mail"   +CC
CC+  command.                                                        +CC
CC+                                                                  +CC
CC+------------------------------------------------------------------+CC
CC+------------------------------------------------------------------+CC
      real function rdl(x,y,z,ier)  
c***begin prologue  rdl
c***date written   790801   (yymmdd)
c***revision date  831228   (yymmdd)
c***category no.  c14
c***keywordls  duplication theorem,elliptic integral,incomplete,complete,
c             integral of the second kind,taylor series
c***author  carlson, b.c., ames laboratory-doe
c             iowa state university, ames, iowa  50011
c           notis, e.m., ames laboratory-doe
c             iowa state university, ames, iowa  50011
c           pexton, r.l., lawrence livermore national laboratory
c             livermore, california  94550
c***purpose  incomplete or complete elliptic integral of the 2nd kind.
c            for x and y nonnegative, x+y and z positive,
c            rdl(x,y,z) = integral from zero to infinity of
c                                -1/2     -1/2     -3/2
c                      (3/2)(t+x)    (t+y)    (t+z)    dt.
c            if x or y is zero, the integral is complete.
c***description
c
c   1.     rdl
c          evaluate an incomplete (or complete) elliptic integral
c          of the second kind
c          standardl fortran function routine
c          single precision version
c          the routine calculates an approximation result to
c          rdl(x,y,z) = integral from zero to infinity of
c                              -1/2     -1/2     -3/2
c                    (3/2)(t+x)    (t+y)    (t+z)    dt,
c          where x and y are nonnegative, x + y is positive, and z is
c          positive.  if x or y is zero, the integral is complete.
c          the duplication theorem is iterated until the variables are
c          nearly equal, and the function is then expanded in taylor
c          series to fifth order.
c
c   2.     calling sequence
c
c          rdl( x, y, z, ier )
c
c          parameters on entry
c          values assigned by the calling routine
c
c          x      - single precision, nonnegative variable
c
c          y      - single precision, nonnegative variable
c
c                   x + y is positive
c
c          z      - real, positive variable
c
c
c
c          on return     (values assigned by the rdl routine)
c
c          rdl     - real approximation to the integral
c
c
c          ier    - integer
c
c                   ier = 0 normal and reliable termination of the
c                           routine.  it is assumed that the requested
c                           accuracy has been achieved.
c
c                   ier >  0 abnormal termination of the routine
c
c
c          x, y, z are unaltered.
c
c   3.    error messages
c
c         value of ier assigned by the rdl routine
c
c                  value assigned         error message printed
c                  ier = 1                amin1(x,y) .lt. 0.0e0
c                      = 2                amin1(x + y, z ) .lt. lolim
c                      = 3                amax1(x,y,z) .gt. uplim
c
c
c   4.     control parameters
c
c                  values of lolim,uplim,and errtol are set by the
c                  routine.
c
c          lolim and uplim determine the valid range of x, y, and z
c
c          lolim  - lower limit of valid arguments
c
c                    not less  than 2 / (machine maximum) ** (2/3).
c
c          uplim  - upper limit of valid arguments
c
c                    not greater than (0.1e0 * errtol / machine
c                    minimum) ** (2/3), where errtol is described below.
c                    in the following table it is assumed that errtol
c                    will never be chosen smaller than 1.0e-5.
c
c
c                    acceptable values for:   lolim      uplim
c                    ibm 360/370 series   :   6.0e-51     1.0e+48
c                    cdc 6000/7000 series :   5.0e-215    2.0e+191
c                    univac 1100 series   :   1.0e-25     2.0e+21
c                    cray                 :   3.0e-1644   1.69e+1640
c                    vax 11 series        :   1.0e-25     4.5e+21
c
c
c          errtol determines the accuracy of the answer
c
c                 the value assigned by the routine will result
c                 in solution precision within 1-2 decimals of
c                 "machine precision".
c
c          errtol    relative error due to truncation is less than
c                    3 * errtol ** 6 / (1-errtol) ** 3/2.
c
c
c
c              the accuracy of the computed approximation to the inte-
c              gral can be controlled by choosing the value of errtol.
c              truncation of a taylor series after terms of fifth ordler
c              introduces an error less than the amount shown in the
c              second column of the following table for each value of
c              errtol in the first column.  in addition to the trunca-
c              tion error there will be round-off error, but in prac-
c              tice the total error from both sources is usually less
c              than the amount given in the table.
c
c
c
c
c          sample choices:  errtol   relative truncation
c                                    error less than
c                           1.0e-3    4.0e-18
c                           3.0e-3    3.0e-15
c                           1.0e-2    4.0e-12
c                           3.0e-2    3.0e-9
c                           1.0e-1    4.0e-6
c
c
c                    decreasing errtol by a factor of 10 yields six more
c                    decimal digits of accuracy at the expense of one or
c                    two more iterations of the duplication theorem.
c***long description
c
c   rdl special comments
c
c
c
c          check: rdl(x,y,z) + rdl(y,z,x) + rdl(z,x,y)
c          = 3 /  sqrt(x * y * z), where x, y, and z are positive.
c
c
c          on input:
c
c          x, y, and z are the variables in the integral rdl(x,y,z).
c
c
c          on output:
c
c
c          x, y, and z are unaltered.
c
c
c
c          ********************************************************
c
c           warning: changes in the program may improve speed at the
c          expense of robustness.
c
c
c
c    -------------------------------------------------------------------
c
c
c   special functions via rdl and rfl
c
c
c                  legendre form of elliptic integral of 2nd kind
c                  ----------------------------------------------
c
c
c                                             2         2   2
c                  e(phi,k) = sin(phi) rfl(cos (phi),1-k sin (phi),1) -
c
c                       2   3             2         2   2
c                  -1/3k sin (phi) rdl(cos (phi),1-k sin (phi),1)
c
c
c                                  2          2         2
c                  e(k) = rfl(0,1-k ,1) - 1/3k rdl(0,1-k ,1)
c
c
c
c                  bulirsch form of elliptic integral of 2nd kind
c                  ----------------------------------------------
c
c                                               2 2    2
c                  el2(x,kc,a,b) = ax rfl(1,1+kc x ,1+x ) +
c
c                                                       2 2    2
c                                 +1/3(b-a) x rdl(1,1+kc x ,1+x )
c
c
c
c                  legendre form of alternative elliptic integral of 2nd
c                  -----------------------------------------------------
c                        kind
c                        ----
c
c                            q     2       2   2  -1/2
c                  d(q,k) = int sin p  (1-k sin p)     dp
c                            0
c
c
c
c                                   3           2     2   2
c                  d(q,k) = 1/3 (sin q)  rdl(cos q,1-k sin q,1)
c
c
c
c
c
c                  lemniscate constant b
c                  ---------------------
c
c
c
c                       1    2    4 -1/2
c                  b = int  s (1-s )    ds
c                       0
c
c
c                  b = 1/3 rdl (0,2,1)
c
c
c
c
c                  heuman's lambda function
c                  ------------------------
c
c
c
c                  pi/2 lambda0(a,b) =
c
c                                        2              2
c                     = sin(b) (rfl(0,cos (a),1)-1/3 sin (a) *
c
c                                2               2         2       2
c                      *rdl(0,cos (a),1)) rfl(cos (b),1-cos (a) sin (b),1)
c
c                             2       3             2
c                     -1/3 cos (a) sin (b) rfl(0,cos (a),1) *
c
c                              2         2       2
c                      *rdl(cos (b),1-cos (a) sin (b),1)
c
c
c
c                  jacobi zeta function
c                  --------------------
c
c
c                             2                 2       2   2
c                  z(b,k) = (k/3) sin(b) rfl(cos (b),1-k sin (b),1)
c
c
c                                       2             2
c                             *rdl(0,1-k ,1)/rfl(0,1-k ,1)
c
c                               2       3           2       2   2
c                            -(k /3) sin (b) rdl(cos (b),1-k sin (b),1)
c
c
c    -------------------------------------------------------------------
c          subroutines or functions needed
c              - errpex
c              - r1pext
c              - fortran abs,amax1,amin1,sqrt
c***references  carlson, b.c. and notis,e .m.
c                 algorithms for incomplete elliptic integrals
c                 acm transactions on mathematical software,vol.7,no.3,
c                 sept, 1981, pages 398-403
c               carlson, b.c.
c                 computing elliptic integrals by duplication
c                 numer. math. 33, (1979), 1-16
c               carlson, b.c.
c                 elliptic integrals of the first kind
c                 siam j. math. anal. 8 (1977), 231-242
c***routines called  r1pext,errpex
c***end prologue  rdl
      integer ier,itodo
      real lolim, uplim, epslon, errtol
      real c1, c2, c3, c4, ea, eb, ec, ed, ef, lamda
      real mu, power4, sigma, s1, s2, x, xn, xndev
      real xnroot, y, yn, yndev, ynroot, z, zn, zndev,
     * znroot
c
c
      character msg*80
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c     the original routine is designed for cft
c
c
c     for  civic  users replace the above character statement
c
c
c        '  character msg*80  '   with the following :
c
c
c           dimension msg(10)
c
c     for  chat  users replace the above character statement
c
c
c        '  character msg*80  '   with the following :
c
c
c           dimension msg(8)
c
c
c      other changes are listed subsequently.
c
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
      data itodo/1/
c
c
c
c***first executable statement  rdl
      if(itodo.eq.1)then
c
c
      errtol=(r1pext(3)/3.0e0)**(1.0e0/6.0e0)
c
c
      lolim = 2.0e0/(r1pext(2))**(2.0e0/3.0e0)
c
      uplim=r1pext(1)**(1.0e0/3.0e0)
      uplim=(0.10e0*errtol)**(1.0e0/3.0e0)/uplim
      uplim=uplim**2.0e0
c
c
      c1 = 3.0e0/14.0e0
      c2 = 1.0e0/6.0e0
      c3 = 9.0e0/22.0e0
      c4 = 3.0e0/26.0e0
c
c
      itodo=0
c
      end if
c
c
c
c         call error handler if necessary.
c
c
c
    5 rdl=0.0e0
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c     for civic users insert the following three statements:
c
c
c          do 9191 igogo=1,10
c
c          msg(igogo)='        '
c
c9191 continue
c
c
c     for chat users insert the following three statements:
c
c
c          do 9191 igogo=1,8
c
c          msg(igogo)='          '
c
c9191 continue
c
c
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
      if( amin1(x,y).lt.0.0e0) then
      ier=1
      msg=  'rdl - error: amin1(x,y).lt.0.0 (where x=r1, y=r2)'
      call errpex (msg,ier,x,y,0.0,0.0,0.0)
      return
      end if
      if (amin1(x+y,z).lt.lolim) then
      ier=2
      msg=  'rdl - error: amin1(x+y,z).lt.lolim (where x=r1, y=r2)'
      call errpex (msg,ier,x,y,0.0,0.0,0.0)
      msg=  'rdl - error: amin1(x+y,z).lt.lolim (where z=r1, lolim=r2)'
      call errpex (msg,ier,z,lolim,0.0,0.0,0.0)
      return
      end if
      if (amax1(x,y,z).gt.uplim) then
      ier=3
      msg=  'rdl - error: amax1(x,y,z).gt.uplim (where x=r1, y=r2)'
      call errpex (msg,ier,x,y,0.0,0.0,0.0)
      msg=  'rdl - error: amax1(x,y,z).gt.uplim (where z=r1, uplim=r2)'
      call errpex (msg,ier,z,lolim,0.0,0.0,0.0)
      return
      end if
c
c
c
   20 ier = 0
      xn = x
      yn = y
      zn = z
      sigma = 0.0e0
      power4 = 1.0e0
c
c
c
   30 mu = (xn+yn+3.0e0*zn)*0.20e0
      xndev = (mu-xn)/mu
      yndev = (mu-yn)/mu
      zndev = (mu-zn)/mu
      epslon = amax1( abs(xndev), abs(yndev), abs(zndev))
      if (epslon.lt.errtol) go to 40
      xnroot =  sqrt(xn)
      ynroot =  sqrt(yn)
      znroot =  sqrt(zn)
      lamda = xnroot*(ynroot+znroot) + ynroot*znroot
      sigma = sigma + power4/(znroot*(zn+lamda))
      power4 = power4*0.250e0
      xn = (xn+lamda)*0.250e0
      yn = (yn+lamda)*0.250e0
      zn = (zn+lamda)*0.250e0
      go to 30
c
c
c
   40 ea = xndev*yndev
      eb = zndev*zndev
      ec = ea - eb
      ed = ea - 6.0e0*eb
      ef = ed + ec + ec
      s1 = ed*(-c1+0.250e0*c3*ed-1.50e0*c4*zndev*ef)
      s2 = zndev*(c2*ef+zndev*(-c3*ec+zndev*c4*ea))
      rdl = 3.0e0*sigma + power4*(1.0e0+s1+s2)/(mu* sqrt(mu))
c
   50 return
c
c
c
      end

      real function rjl(x,y,z,p,ier)       
c***begin prologue  rjl
c***date written   790801   (yymmdd)
c***revision date  831229   (yymmdd)
c***category no.  c14
c***keywords  duplication theorem,elliptic integral,incomplete,complete,
c             integral of the third kind,taylor series
c***author  carlson, b.c., ames laboratory-doe
c             iowa state university, ames, iowa  50011
c           notis, e.m., ames laboratory-doe
c             iowa state university, ames, iowa  50011
c           pexton, r.l., lawrence livermore national laboratory
c             livermore, california  94550
c***purpose  incomplete or complete elliptic integral of the 3rd kind.
c            for x, y, and z nonnegative, at most one of them zero, and
c            p positive, rjl(x,y,z,p) = integral from zero to infinity of
c                                  -1/2     -1/2     -1/2     -1
c                        (3/2)(t+x)    (t+y)    (t+z)    (t+p)  dt,
c            if x or y or z is 0, the integral is complete.
c***description
c
c   1.     rjl
c          standard fortran function routine
c          single precision version
c          the routine calculates an approximation result to
c          rjl(x,y,z,p) = integral from zero to infinity of
c
c                                -1/2     -1/2     -1/2     -1
c                      (3/2)(t+x)    (t+y)    (t+z)    (t+p)  dt,
c
c          where x, y, and z are nonnegative, at most one of them is
c          zero, and p is positive.  if x or y or z is zero, the
c          integral is complete.  the duplication theorem is iterated
c          until the variables are nearly equal, and the function is
c          then expanded in taylor series to fifth order.
c
c
c   2.     calling sequence
c          rjl( x, y, z, p, ier )
c
c          parameters on entry
c          values assigned by the calling routine
c
c          x      - single precision, nonnegative variable
c
c          y      - single precision, nonnegative variable
c
c          z      - single precision, nonnegative variable
c
c          p      - single precision, positive variable
c
c
c          on  return     (values assigned by the rjl routine)
c
c          rjl     - single precision approximation to the integral
c
c          ier    - integer
c
c                   ier = 0 normal and reliable termination of the
c                           routine.  it is assumed that the requested
c                           accuracy has been achieved.
c
c                   ier >  0 abnormal termination of the routine
c
c
c          x, y, z, p are unaltered.
c
c
c   3.    error messages
c
c         value of ier assigned by the rjl routine
c
c                  value assigned        error message printed
c                  ier = 1               amin1(x,y,z) .lt. 0.0e0
c                      = 2               amin1(x+y,x+z,y+z,p) .lt. lolim
c                      = 3               amax1(x,y,z,p) .gt. uplim
c
c
c
c   4.     control parameters
c
c                  values of lolim,uplim,and errtol are set by the
c                  routine.
c
c
c          lolim and uplim determine the valid range of x y, z, and p
c
c          lolim is not less than the cube root of the value
c          of lolim used in the routine for rcl.
c
c          uplim is not greater than 0.3 times the cube root of
c          the value of uplim used in the routine for rcl.
c
c
c                     acceptable values for:   lolim      uplim
c                     ibm 360/370 series   :   2.0e-26     3.0e+24
c                     cdc 6000/7000 series :   5.0e-98     3.0e+106
c                     univac 1100 series   :   5.0e-13     6.0e+11
c                     cray                 :   1.32e-822   1.4e+821
c                     vax 11 series        :   2.5e-13     9.0e+11
c
c
c
c          errtol determines the accuracy of the answer
c
c                 the value assigned by the routine will result
c                 in solution precision within 1-2 decimals of
c                 "machine precision".
c
c
c
c
c          relative error due to truncation of the series for rjl
c          is less than 3 * errtol ** 6 / (1 - errtol) ** 3/2.
c
c
c
c              the accuracy of the computed approximation to the inte-
c              gral can be controlled by choosing the value of errtol.
c              truncation of a taylor series after terms of fifth order
c              introduces an error less than the amount shown in the
c              second column of the following table for each value of
c              errtol in the first column.  in addition to the trunca-
c              tion error there will be round-off error, but in prac-
c              tice the total error from both sources is usually less
c              than the amount given in the table.
c
c
c
c          sample choices:  errtol   relative truncation
c                                    error less than
c                           1.0e-3    4.0e-18
c                           3.0e-3    3.0e-15
c                           1.0e-2    4.0e-12
c                           3.0e-2    3.0e-9
c                           1.0e-1    4.0e-6
c
c                    decreasing errtol by a factor of 10 yields six more
c                    decimal digits of accuracy at the expense of one or
c                    two more iterations of the duplication theorem.
c***long description
c
c   rjl special comments
c
c
c          check by addition theorem: rjl(x,x+z,x+w,x+p)
c          + rjl(y,y+z,y+w,y+p) + (a-b) * rjl(a,b,b,a) + 3 / sqrt(a)
c          = rjl(0,z,w,p), where x,y,z,w,p are positive and x * y
c          = z * w,  a = p * p * (x+y+z+w),  b = p * (p+x) * (p+y),
c          and b - a = p * (p-z) * (p-w).  the sum of the third and
c          fourth terms on the left side is 3 * rcl(a,b).
c
c
c          on input:
c
c          x, y, z, and p are the variables in the integral rjl(x,y,z,p).
c
c
c          on output:
c
c
c          x, y, z, and p are unaltered.
c
c          ********************************************************
c
c          warning: changes in the program may improve speed at the
c          expense of robustness.
c
c ------------------------------------------------------------
c
c
c   special functions via rjl,rdl, and rfl
c
c
c                  legendre form of elliptic integral of 3rd kind
c                  ----------------------------------------------
c
c
c                               phi         2         -1
c                  p(phi,k,n) = int (1+n sin (theta) )   *
c                                0
c
c                                      2    2         -1/2
c                                 *(1-k  sin (theta) )     d theta
c
c
c                                               2          2   2
c                            = sin (phi) rfl(cos (phi), 1-k sin (phi),1)
c
c                                     3             2         2   2
c                             -n/3 sin (phi) rjl(cos (phi),1-k sin (phi),
c
c                                      2
c                             1,1+n sin (phi))
c
c
c
c                  bulirsch form of elliptic integral of 3rd kind
c                  ----------------------------------------------
c
c
c                                            2 2    2
c                  el3(x,kc,p) = x rfl(1,1+kc x ,1+x ) +
c
c                                          3           2 2    2     2
c                               +1/3(1-p) x  rjl(1,1+kc x ,1+x ,1+px )
c
c
c                                            2
c                  cel(kc,p,a,b) = a rfl(0,kc ,1) +
c
c
c                                 +1/3(b-pa) rjl(0,kc ,1,p)
c
c
c
c
c                  heuman's lambda function
c                  ------------------------
c
c
c                                 2                     2      2    1/2
c                  l(a,b,p) = (cos(a)sin(b)cos(b)/(1-cos (a)sin (b))   )
c
c                                            2         2       2
c                            *(sin(p) rfl(cos (p),1-sin (a) sin (p),1)
c
c                                 2       3            2       2
c                            +(sin (a) sin (p)/(3(1-cos (a) sin (b))
c
c                                    2         2       2
c                            *rjl(cos (p),1-sin (a) sin (p),1,1-
c
c                                2       2          2       2
c                            -sin (a) sin (p)/(1-cos (a) sin (b))))
c
c
c
c
c                  pi/2 lambda0(a,b) =l(a,b,pi/2) =
c
c                                        2              2
c                     = sin(b) (rfl(0,cos (a),1)-1/3 sin (a) *
c
c                               2               2         2       2
c                     *rdl(0,cos (a),1)) rfl(cos (b),1-cos (a) sin (b),1)
c
c                             2       3             2
c                     -1/3 cos (a) sin (b) rfl(0,cos (a),1) *
c
c                             2         2       2
c                     *rdl(cos (b),1-cos (a) sin (b),1)
c
c
c
c                  jacobi zeta function
c                  --------------------
c
c
c                             2                     2   2    1/2
c                  z(b,k) = (k/3) sin(b) cos(b) (1-k sin (b))
c
c
c                                       2      2   2                 2
c                             *rjl(0,1-k ,1,1-k sin (b)) / rfl (0,1-k ,1)
c
c
c    -------------------------------------------------------------------
c
c          subroutine or functions needed
c              - rcl
c              - errpex
c              - r1pext
c              - fortran abs,amax1,amin1,sqrt
c***references  carlson, b.c. and notis, e.m.
c                 algorithms for incomplete elliptic integrals
c                 acm transactions on mathematical software,vol.7,no.3,
c                 sept, 1981, pages 398-403
c               carlson, b.c.
c                 computing elliptic integrals by duplication
c                 numer. math. 33, (1979), 1-16
c               carlson, b.c.
c                 elliptic integrals of the first kind
c                 siam j. math. anal. 8 (1977), 231-242
c***routines called  r1pext,rcl,errpex
c***end prologue  rjl
      integer ier,itodo
      real alfa, beta, c1, c2, c3, c4, ea, eb, ec, e2, e3
      real lolim, uplim, epslon, errtol, etolrcl
      real lamda, mu, p, pn, pndev
      real power4, rcl, sigma, s1, s2, s3, x, xn, xndev
      real xnroot, y, yn, yndev, ynroot, z, zn, zndev,
     * znroot
c
c
      character msg*80
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c     the original routine is designed for cft
c
c
c     for  civic  users replace the above character statement
c
c
c        '  character msg*80  '   with the following :
c
c
c           dimension msg(10)
c
c
c     for  chat  users replace the above character statement
c
c
c        '  character msg*80  '   with the following :
c
c
c           dimension msg(8)
c
c
c      other changes are listed subsequently.
c
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c
c
      data itodo/1/
c
c
c
c***first executable statement  rjl
      if(itodo.eq.1)then
c
c
      errtol=(r1pext(3)/3.0e0)**(1.0e0/6.0e0)
c
c
      lolim =( 5.0e0 * r1pext(1))**(1.0e0/3.0e0)
c
      uplim = 0.30e0*( r1pext(2) / 5.0e0)**(1.0e0/3.0e0)
c
c
c
      c1 = 3.0e0/14.0e0
      c2 = 1.0e0/3.0e0
      c3 = 3.0e0/22.0e0
      c4 = 3.0e0/26.0e0
c
c
c
      itodo=0
c
      end if
c
c
c
c         call error handler if necessary.
c
c
    5 rjl=0.0e0
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c     for civic users insert the following three statements:
c
c
c          do 9191 igogo=1,10
c
c          msg(igogo)='        '
c
c9191 continue
c
c
c     for chat users insert the following three statements:
c
c
c          do 9191 igogo=1,8
c
c          msg(igogo)='          '
c
c9191 continue
c
c
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
      if( amin1(x,y,z).lt.0.0e0) then
      ier=1
      msg='rjl-error: amin1(x,y,z).lt.0.0e0 (where x=r1, y=r2, z=r3)'
      call errpex (msg,ier,x,y,z,p,uplim)
      return
      end if
      if (amin1(x+y,x+z,y+z,p).lt.lolim) then
      ier=2
      msg= 'rjl-error: amin1(x+y,x+z,y+z,p).lt.lolim (where x=r1, y=r2,
     *z=r3, p=r4, lolim=r5)'
      call errpex (msg,ier,x,y,z,p,lolim)
      return
      end if
      if (amax1(x,y,z,p).gt.uplim) then
      ier=3
      msg= 'rjl-error: amax1(x,y,z,p).gt.uplim (where x=r1, y=r2,
     *z=r3, p=r4, uplim=r5)'
      call errpex (msg,ier,x,y,z,p,uplim)
      return
      end if
c
c
c
   20 ier = 0
      xn = x
      yn = y
      zn = z
      pn = p
      sigma = 0.0e0
      power4 = 1.0e0
c
c
c
   30 mu = (xn+yn+zn+pn+pn)*0.20e0
      xndev = (mu-xn)/mu
      yndev = (mu-yn)/mu
      zndev = (mu-zn)/mu
      pndev = (mu-pn)/mu
      epslon = amax1( abs(xndev), abs(yndev), abs(zndev), abs(pndev))
      if (epslon.lt.errtol) go to 40
      xnroot =  sqrt(xn)
      ynroot =  sqrt(yn)
      znroot =  sqrt(zn)
      lamda = xnroot*(ynroot+znroot) + ynroot*znroot
      alfa = pn*(xnroot+ynroot+znroot) + xnroot*ynroot*znroot
      alfa = alfa*alfa
      beta = pn*(pn+lamda)*(pn+lamda)
      sigma = sigma + power4*rcl(alfa,beta,ier)
      power4 = power4*0.250e0
      xn = (xn+lamda)*0.250e0
      yn = (yn+lamda)*0.250e0
      zn = (zn+lamda)*0.250e0
      pn = (pn+lamda)*0.250e0
      go to 30
c
c
c
   40 ea = xndev*(yndev+zndev) + yndev*zndev
      eb = xndev*yndev*zndev
      ec = pndev*pndev
      e2 = ea - 3.0e0*ec
      e3 = eb + 2.0e0*pndev*(ea-ec)
      s1 = 1.0e0 + e2*(-c1+0.750e0*c3*e2-1.50e0*c4*e3)
      s2 = eb*(0.50e0*c2+pndev*(-c3-c3+pndev*c4))
      s3 = pndev*ea*(c2-pndev*c3) - c2*pndev*ec
      rjl = 3.0e0*sigma + power4*(s1+s2+s3)/(mu* sqrt(mu))
c
c
   50 return
c
c
c
      end

      real function rcl(x,y,ier)    
      integer ier,itodo
      real c1, c2, errtol, lamda, lolim
      real mu, s, sn, uplim, x, xn, y, yn
c
c***begin prologue  rcl
c***date written   790801   (yymmdd)
c***revision date  830622   (yymmdd)
c***category no.  c14
c***keywords  duplication theorem,elementary functions,
c             elliptic integrals,taylor series
c
c
c***author  carlson, b.c., ames laboratory-doe
c             iowa state university, ames, iowa  50011
c           notis, e.m., ames laboratory-doe,
c             iowa state university, ames, iowa  50011
c           pexton, r.l., lawrence livermore national laboratory
c             livermore, california  94550
c***purpose  the routine calculates an approximation result to
c            rcl(x,y) = integral from zero to infinity of
c                              -1/2     -1
c                    (1/2)(t+x)    (t+y)  dt,
c            where x is nonnegative and y is positive.
c***description
c
c   1.     rcl
c          standard fortran function routine
c          single precision version
c          the routine calculates an approximation result to
c          rcl(x,y) = integral from zero to infinity of
c
c                              -1/2     -1
c                    (1/2)(t+x)    (t+y)  dt,
c
c          where x is nonnegative and y is positive.  the duplication
c          theorem is iterated until the variables are nearly equal,
c          and the function is then expanded in taylor series to fifth
c          order.  logarithmic, inverse circular, and inverse hyper-
c          bolic functions can be expressed in terms of rcl.
c
c
c   2.     calling sequence
c          rcl( x, y, ier )
c
c          parameters on entry
c          values assigned by the calling routine
c
c          x      - single precision, nonnegative variable
c
c          y      - single precision, positive variable
c
c
c
c          on return  (values assigned by the rcl routine)
c
c          rcl     - single precision approximation to the integral
c
c          ier    - integer to indicate normal or abnormal termination.
c
c                     ier = 0 normal and reliable termination of the
c                             routine.  it is assumed that the requested
c                             accuracy has been achieved.
c
c                     ier > 0 abnormal termination of the routine
c
c          x and y are unaltered.
c
c
c   3.    error messages
c
c         value of ier assigned by the rcl routine
c
c                  value assigned         error message printed
c                  ier = 1                x.lt.0.0e0.or.y.le.0.0e0
c                      = 2                x+y.lt.lolim
c                      = 3                amax1(x,y) .gt. uplim
c
c
c   4.     control parameters
c
c                  values of lolim,uplim,and errtol are set by the
c                  routine.
c
c          lolim and uplim determine the valid range of x and y
c
c          lolim  - lower limit of valid arguments
c
c                   not less  than 5 * (machine minimum)  .
c
c          uplim  - upper limit of valid arguments
c
c                   not greater than (machine maximum) / 5 .
c
c
c                     acceptable values for:   lolim       uplim
c                     ibm 360/370 series   :   3.0e-78     1.0e+75
c                     cdc 6000/7000 series :   1.0e-292    1.0e+321
c                     univac 1100 series   :   1.0e-37     1.0e+37
c                     cray                 :   2.3e-2466   1.09e+2465
c                     vax 11 series        :   1.5e-38     3.0e+37
c
c          errtol determines the accuracy of the answer
c
c                 the value assigned by the routine will result
c                 in solution precision within 1-2 decimals of
c                 "machine precision".
c
c
c          errtol  - relative error due to truncation is less than
c                    16 * errtol ** 6 / (1 - 2 * errtol).
c
c
c              the accuracy of the computed approximation to the inte-
c              gral can be controlled by choosing the value of errtol.
c              truncation of a taylor series after terms of fifth order
c              introduces an error less than the amount shown in the
c              second column of the following table for each value of
c              errtol in the first column.  in addition to the trunca-
c              tion error there will be round-off error, but in prac-
c              tice the total error from both sources is usually less
c              than the amount given in the table.
c
c
c
c          sample choices:  errtol   relative truncation
c                                    error less than
c                           1.0e-3    2.0e-17
c                           3.0e-3    2.0e-14
c                           1.0e-2    2.0e-11
c                           3.0e-2    2.0e-8
c                           1.0e-1    2.0e-5
c
c
c                    decreasing errtol by a factor of 10 yields six more
c                    decimal digits of accuracy at the expense of one or
c                    two more iterations of the duplication theorem.
c***long description
c
c   rcl special comments
c
c
c
c
c                  check: rcl(x,x+z) + rcl(y,y+z) = rcl(0,z)
c
c                  where x, y, and z are positive and x * y = z * z
c
c
c          on input:
c
c          x and y are the variables in the integral rcl(x,y).
c
c          on output:
c
c          x and y are unaltered.
c
c
c
c                    rcl(0,1/4)=rcl(1/16,1/8)=pi=3.14159...
c
c                    rcl(9/4,2)=ln(2)
c
c
c
c          ********************************************************
c
c          warning: changes in the program may improve speed at the
c          expense of robustness.
c
c
c   --------------------------------------------------------------------
c
c   special functions via rcl
c
c
c
c                  ln x                x .gt. 0
c
c                                             2
c                  ln(x) = (x-1) rcl(((1+x)/2)  , x )
c
c
c   --------------------------------------------------------------------
c
c                   arcsin x            -1 .le. x .le. 1
c
c                                        2
c                   arcsin x = x rcl (1-x  ,1 )
c
c   --------------------------------------------------------------------
c
c                   arccos x            0 .le. x .le. 1
c
c
c                                      2        2
c                   arccos x = sqrt(1-x ) rcl(x  ,1 )
c
c   --------------------------------------------------------------------
c
c                   arctan x            -inf .lt. x .lt. +inf
c
c                                         2
c                   arctan x = x rcl(1,1+x  )
c
c   --------------------------------------------------------------------
c
c                   arccot x            0 .le. x .lt. inf
c
c                                   2   2
c                   arccot x = rcl(x  ,x +1 )
c
c   --------------------------------------------------------------------
c
c                   arcsinh x           -inf .lt. x .lt. +inf
c
c                                        2
c                   arcsinh x = x rcl(1+x  ,1 )
c
c   --------------------------------------------------------------------
c
c                   arccosh x           x .ge. 1
c
c                                     2         2
c                   arccosh x = sqrt(x -1) rcl(x  ,1 )
c
c   --------------------------------------------------------------------
c
c                   arctanh x           -1 .lt. x .lt. 1
c
c                                          2
c                   arctanh x = x rcl(1,1-x  )
c
c   --------------------------------------------------------------------
c
c                   arccoth x           x .gt. 1
c
c                                    2   2
c                   arccoth x = rcl(x  ,x -1 )
c
c   --------------------------------------------------------------------
c
c          subroutines or functions needed
c              - errpex
c              - r1pext
c              - fortran abs,amax1,sqrt
c***references  carlson, b.c. and notis, e.m.
c                 algorithms for incomplete elliptic integrals
c                 acm transactions on mathematical software,vol.7,no.3,
c                 sept, 1981, pages 398-403
c               carlson, b.c.
c                 computing elliptic integrals by duplication
c                 numer. math. 33, (1979), 1-16
c               carlson, b.c.
c                 elliptic integrals of the first kind
c                 siam j. math. anal. 8 (1977), 231-242
c***routines called  r1pext,errpex
c
      character msg*80
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c     the original routine is designed for cft
c
c
c     for  civic  users replace the above character statement
c
c
c        '  character msg*80  '   with the following :
c
c
c           dimension msg(10)
c
c
c     for  chat  users replace the above character statement
c
c
c        '  character msg*80  '   with the following :
c
c
c           dimension msg(8)
c
c
c      other changes are listed subsequently.
c
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
      data itodo/1/
c***first executable statement  rcl
      if(itodo.eq.1)then
c
c
      errtol=(r1pext(3)/16.0)**(1.0/6.0)
c
c
      lolim = 5.0e0 * r1pext(1)
      uplim = r1pext(2) / 5.0e0
c
c
      c1 = 1.0e0/7.0e0
      c2 = 9.0e0/22.0e0
c
c
c
      itodo=0
c
      end if
c
c
c         call error handler if necessary.
c
c
    5 rcl=0.0e0
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
c     for civic users insert the following three statements:
c
c
c          do 9191 igogo=1,10
c
c          msg(igogo)='        '
c
c9191 continue
c
c
c     for chat users insert the following three statements:
c
c
c          do 9191 igogo=1,8
c
c          msg(igogo)='          '
c
c9191 continue
c
c
c
c
c###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###   ###
c
c
      if (x.lt.0.0e0.or.y.le.0.0e0) then
      ier=1
      msg='rcl - error: x.lt.0.0e0.or.y.le.0.0e0 (where x=r1, y=r2)'
      call errpex (msg,ier,x,y,0.0,0.0,0.0)
      return
      end if
      if (x+y.lt.lolim)   then
      ier=2
      msg='rcl - error: x+y.lt.lolim (where x=r1, y=r2)'
      call errpex (msg,ier,x,y,0.0,0.0,0.0)
      msg='rcl - error: x+y.lt.lolim (where lolim=r1)'
      call errpex (msg,ier,lolim,0.0,0.0,0.0,0.0)
      return
      end if
      if (amax1(x,y).gt.uplim) then
      ier=3
      msg='rcl - error: amax1(x,y).gt.uplim (where x=r1, y=r2)'
      call errpex (msg,ier,x,y,0.0,0.0,0.0)
      msg='rcl - error: amax1(x,y).gt.uplim (where uplim=r1)'
      call errpex (msg,ier,uplim,0.0,0.0,0.0,0.0)
      return
      end if
c
c
c
   20 ier = 0
      xn = x
      yn = y
c
c
c
   30 mu = (xn+yn+yn)/3.0e0
      sn = (yn+mu)/mu - 2.0e0
      if ( abs(sn).lt.errtol) go to 40
      lamda = 2.0e0* sqrt(xn)* sqrt(yn) + yn
      xn = (xn+lamda)*0.250e0
      yn = (yn+lamda)*0.250e0
      go to 30
c
   40 s = sn*sn*(0.30e0+sn*(c1+sn*(0.3750e0+sn*c2)))
      rcl = (1.0e0+s)/ sqrt(mu)
c
c
c
   50 return
c
c
c
c
      end