!***********************************************************************************************************************************
!
!                                                               A G M
!
!  Program:      AGM
!
!  Programmer:   David G. Simpson
!                NASA Goddard Space Flight Center
!                Greenbelt, Maryland  20771
!
!  Date:         November 22, 2025
!
!  Language:     Fortran-90
!
!  Version:      1.00a
!
!  Description:  Function AGM() to compute arithmetic-geometric mean, along with a main test program.
!
!***********************************************************************************************************************************

      PROGRAM AGM_TEST

      IMPLICIT NONE

      DOUBLE PRECISION :: AGM


      PRINT *, ' agm(24,6)         = ', AGM(24.0D0, 6.0D0, 1.0D-12)       ! should be 13.458171...
      PRINT *, ' Gauss''s constant  =  ', 1.0D0/(AGM(1.0D0, SQRT(2.0D0)))  ! should be  0.83462684167...

      STOP

      END PROGRAM AGM_TEST





!***********************************************************************************************************************************
!  AGM
!
!  Returns the arithmetic-geometric mean of the two arguments A and B.
!
!  The optional third argument TOL, if specified, is the convergence tolerance.  If not specified, TOL defaults to 1.0D-12.
!  The returned value is the arithmetic-geometric mean, agm(A,B).
!
!***********************************************************************************************************************************

      RECURSIVE FUNCTION AGM (A, B, TOL) RESULT (ANS)

      IMPLICIT NONE

      DOUBLE PRECISION, INTENT(IN) :: A, B
      DOUBLE PRECISION, INTENT(IN), OPTIONAL :: TOL
      DOUBLE PRECISION :: ANS, ARITH_MEAN, GEOM_MEAN, CUTOFF

      IF (PRESENT(TOL)) THEN
         CUTOFF = TOL
      ELSE
         CUTOFF = 1.0D-12
      END IF

      ARITH_MEAN = 0.5D0*(A+B)
      GEOM_MEAN = SQRT(A*B)
      IF (ABS(ARITH_MEAN-GEOM_MEAN) .LT. CUTOFF) THEN
         ANS = A
      ELSE
         ANS = AGM(ARITH_MEAN, GEOM_MEAN, CUTOFF)
      END IF

      RETURN

      END FUNCTION AGM