*     G01MUF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          MODE
      PARAMETER        (MODE=0)
*     .. Local Scalars ..
      DOUBLE PRECISION BETA2, C1, C2, RKAPPA, X, XL, XU, Y
      INTEGER          IFAIL
*     .. External Functions ..
      DOUBLE PRECISION G01MUF, X02ALF
      EXTERNAL         G01MUF, X02ALF
*     .. Local Arrays ..
      DOUBLE PRECISION WORK(322)
*     .. External Subroutines ..
      EXTERNAL         G01ZUF
*     .. Executable Statements ..
      WRITE (NOUT,*) 'G01MUF Example Program Results'
*     Skip heading in data file
      READ (NIN,*)
      C1 = -X02ALF()
      C2 = -X02ALF()
      WRITE (NOUT,*)
      WRITE (NOUT,*) '  X     RKAPPA    BETA2        Y        IFAIL'
      WRITE (NOUT,*)
   20 READ (NIN,*,END=40) X, RKAPPA, BETA2
      IF ((RKAPPA.NE.C1) .OR. (BETA2.NE.C2)) THEN
*
*        Initialize array WORK before the the first call to G01MUF and
*        on subsequent calls when RKAPPA or BETA2 has changed.
*
         IFAIL = 1
*
         CALL G01ZUF(RKAPPA,BETA2,MODE,XL,XU,WORK,IFAIL)
*
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99999) ' ** G01ZUF returned with IFAIL = ',
     +        IFAIL
            GO TO 40
         END IF
*
      END IF
*
*     Compute the value of the Vavilov density function
*
      IFAIL = 1
*
      Y = G01MUF(X,WORK,IFAIL)
*
      IF (IFAIL.EQ.0) THEN
         WRITE (NOUT,99998) X, RKAPPA, BETA2, Y, IFAIL
         C1 = RKAPPA
         C2 = BETA2
         GO TO 20
      ELSE
         WRITE (NOUT,99999) ' ** G01MUF returned with IFAIL = ', IFAIL
      END IF
*
   40 CONTINUE
*
99999 FORMAT (1X,A,I5)
99998 FORMAT (1X,F4.1,5X,F4.1,5X,F4.1,3X,1P,E12.4,I6)
      END