*     D06BAF Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NEDMX, NVMAX, NLINESX, SDCRUS, NCOMPX, NVIMX,
     +                 MAXCAN, LRWORK, LIWORK
      PARAMETER        (NEDMX=1000,NVMAX=5000,NLINESX=50,SDCRUS=1000,
     +                 NCOMPX=5,NVIMX=20,MAXCAN=10000,
     +                 LRWORK=12*NVMAX+3*MAXCAN+15,
     +                 LIWORK=8*NEDMX+55*NVMAX+MAXCAN+78)
*     .. Local Scalars ..
      DOUBLE PRECISION X0, XA, XB, XMAX, XMIN, Y0, YMAX, YMIN
      INTEGER          I, IFAIL, ITRACE, J, K, NCOMP, NEDGE, NELT,
     +                 NLINES, NPROPA, NV, NVB, NVINT, REFTK
      CHARACTER        PMESH
*     .. Local Arrays ..
      DOUBLE PRECISION COOR(2,NVMAX), COORCH(2,NLINESX), CRUS(2,SDCRUS),
     +                 RATE(NLINESX), RUSER(4), RWORK(LRWORK),
     +                 WEIGHT(NVIMX)
      INTEGER          CONN(3,2*NVMAX+5), EDGE(3,NEDMX), IUSER(1),
     +                 IWORK(LIWORK), LCOMP(NLINESX), LINED(4,NLINESX),
     +                 NLCOMP(NCOMPX)
*     .. External Functions ..
      DOUBLE PRECISION FBND
      EXTERNAL         FBND
*     .. External Subroutines ..
      EXTERNAL         D06ABF, D06ACF, D06BAF
*     .. Intrinsic Functions ..
      INTRINSIC        ABS
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D06BAF Example Program Results'
      WRITE (NOUT,*)
*
*     Skip heading in data file
*
      READ (NIN,*)
*
*     Initialise boundary mesh inputs:
*     the number of line and of the characteristic points of
*     the boundary mesh
*
      READ (NIN,*) NLINES
*
*     The ellipse boundary which envelops the NAg Logo
*     the N, the A and the G
*
      READ (NIN,*) (COORCH(1,J),J=1,NLINES)
      READ (NIN,*) (COORCH(2,J),J=1,NLINES)
*
      READ (NIN,*) (CRUS(1,J),J=1,4)
      READ (NIN,*) (CRUS(2,J),J=1,4)
*
*     The Lines of the boundary mesh
*
      READ (NIN,*) ((LINED(I,J),I=1,4),RATE(J),J=1,NLINES)
*
*     The number of connected components to the boundary
*     and their informations
*
      READ (NIN,*) NCOMP
      J = 1
      DO 20 I = 1, NCOMP
         READ (NIN,*) NLCOMP(I)
*
         READ (NIN,*) (LCOMP(K),K=J,J+ABS(NLCOMP(I))-1)
         J = J + ABS(NLCOMP(I))
   20 CONTINUE
*
      READ (NIN,*) PMESH
*
*     Data passed to the user-supplied function
*
      XMIN = COORCH(1,4)
      XMAX = COORCH(1,2)
      YMIN = COORCH(2,1)
      YMAX = COORCH(2,3)
*
      XA = (XMAX-XMIN)/2.D0
      XB = (YMAX-YMIN)/2.D0
*
      X0 = (XMIN+XMAX)/2.D0
      Y0 = (YMIN+YMAX)/2.D0
*
      RUSER(1) = XA
      RUSER(2) = XB
      RUSER(3) = X0
      RUSER(4) = Y0
      IUSER(1) = 0
*
      ITRACE = -1
*
*     Call to the boundary mesh generator
*
      IFAIL = 1
*
      CALL D06BAF(NLINES,COORCH,LINED,FBND,CRUS,SDCRUS,RATE,NCOMP,
     +            NLCOMP,LCOMP,NVMAX,NEDMX,NVB,COOR,NEDGE,EDGE,ITRACE,
     +            RUSER,IUSER,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
*
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,*)
         WRITE (NOUT,99993) ' ** D06BAF returned with IFAIL = ', IFAIL
         GO TO 160
      ELSE IF (PMESH.EQ.'N') THEN
         WRITE (NOUT,*) 'Boundary mesh characteristics'
         WRITE (NOUT,99999) 'NVB   =', NVB
         WRITE (NOUT,99999) 'NEDGE =', NEDGE
      ELSE IF (PMESH.EQ.'Y') THEN
*
*        Output the mesh to view it using the NAG Graphics Library
*
         WRITE (NOUT,99998) NVB, NEDGE
*
         DO 40 I = 1, NVB
            WRITE (NOUT,99997) I, COOR(1,I), COOR(2,I)
   40    CONTINUE
*
         DO 60 I = 1, NEDGE
            WRITE (NOUT,99996) I, EDGE(1,I), EDGE(2,I), EDGE(3,I)
   60    CONTINUE
      ELSE
         WRITE (NOUT,*) 'Problem with the printing option Y or N'
         GO TO 160
      END IF
*
*     Initialise mesh control parameters
*
      ITRACE = 0
      NPROPA = 1
      NVINT = 0
      IFAIL = 1
*
*     Call to the 2D Delaunay-Voronoi mesh generator
*
      CALL D06ABF(NVB,NVINT,NVMAX,NEDGE,EDGE,NV,NELT,COOR,CONN,WEIGHT,
     +            NPROPA,ITRACE,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
*
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,*)
         WRITE (NOUT,99993) ' ** D06ABF returned with IFAIL = ', IFAIL
         GO TO 160
      ELSE IF (PMESH.EQ.'N') THEN
         WRITE (NOUT,*) 'Complete mesh (via the 2D Delaunay-Voronoi'
         WRITE (NOUT,*) 'mesh generator) characteristics'
         WRITE (NOUT,99999) 'NV   =', NV
         WRITE (NOUT,99999) 'NELT =', NELT
      ELSE IF (PMESH.EQ.'Y') THEN
*
*        Output the mesh to view it using the NAG Graphics Library
*
         WRITE (NOUT,99998) NV, NELT
         DO 80 I = 1, NV
            WRITE (NOUT,99995) COOR(1,I), COOR(2,I)
   80    CONTINUE
*
         REFTK = 0
         DO 100 K = 1, NELT
            WRITE (NOUT,99994) CONN(1,K), CONN(2,K), CONN(3,K), REFTK
  100    CONTINUE
      END IF
*
*     Call to the 2D Advancing front mesh generator
*
      IFAIL = 1
*
      CALL D06ACF(NVB,NVINT,NVMAX,NEDGE,EDGE,NV,NELT,COOR,CONN,WEIGHT,
     +            ITRACE,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
*
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,*)
         WRITE (NOUT,99993) ' ** D06ACF returned with IFAIL = ', IFAIL
         GO TO 160
      ELSE IF (PMESH.EQ.'N') THEN
         WRITE (NOUT,*) 'Complete mesh (via the 2D Advancing front mesh'
         WRITE (NOUT,*) 'generator) characteristics'
         WRITE (NOUT,99999) 'NV   =', NV
         WRITE (NOUT,99999) 'NELT =', NELT
      ELSE IF (PMESH.EQ.'Y') THEN
*
*        Output the mesh to view it using the NAG Graphics Library
*
         WRITE (NOUT,99998) NV, NELT
         DO 120 I = 1, NV
            WRITE (NOUT,99995) COOR(1,I), COOR(2,I)
  120    CONTINUE
*
         REFTK = 0
         DO 140 K = 1, NELT
            WRITE (NOUT,99994) CONN(1,K), CONN(2,K), CONN(3,K), REFTK
  140    CONTINUE
      END IF
*
  160 CONTINUE
*
99999 FORMAT (1X,A,I6)
99998 FORMAT (1X,2I10)
99997 FORMAT (2X,I4,2(2X,E12.6))
99996 FORMAT (1X,4I4)
99995 FORMAT (2(2X,E12.6))
99994 FORMAT (1X,4I10)
99993 FORMAT (1X,A,I5)
      END
*
      DOUBLE PRECISION FUNCTION FBND(I,X,Y,RUSER,IUSER)
*     .. Scalar Arguments ..
      DOUBLE PRECISION X, Y
      INTEGER          I
*     .. Array Arguments ..
      DOUBLE PRECISION RUSER(*)
      INTEGER          IUSER(*)
*     .. Local Scalars ..
      DOUBLE PRECISION RADIUS2, X0, XA, XB, Y0
*     .. Executable Statements ..
      XA = RUSER(1)
      XB = RUSER(2)
      X0 = RUSER(3)
      Y0 = RUSER(4)
*
      FBND = 0.D0
      IF (I.EQ.1) THEN
*
*        line 1,2,3, and 4: ellipse centred in (X0,Y0) with
*        XA and XB as coefficients
*
         FBND = ((X-X0)/XA)**2 + ((Y-Y0)/XB)**2 - 1.D0
      ELSE IF (I.EQ.2) THEN
*
*        line 7, lower arc on letter n, is a circle centred in (X0,Y0)
*        with radius SQRT(RADIUS2)
*
         X0 = 0.5D0
         Y0 = 6.25D0
         RADIUS2 = 20.3125D0
         FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2
      ELSE IF (I.EQ.3) THEN
*
*        line 11, upper arc on letter n, is a circle centred in (X0,Y0)
*        with radius SQRT(RADIUS2)
*
         X0 = 1.0D0
         Y0 = 4.0D0
         RADIUS2 = 9.0D0 + (11.0D0-Y0)**2
         FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2
      ELSE IF (I.EQ.4) THEN
*
*        line 15, upper arc on letter a, is a circle centred in (X0,Y0)
*        with radius SQRT(RADIUS2) touching point (5,11).
*
         X0 = 8.5D0
         Y0 = 2.75D0
         RADIUS2 = (X0-5.0D0)**2 + (11.0D0-Y0)**2
         FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2
      ELSE IF (I.EQ.5) THEN
*
*        line 25, lower arc on hat of 'a', is a circle centred in (X0,Y0)
*        with radius SQRT(RADIUS2) touching point (11,10).
*
         X0 = 8.5D0
         Y0 = 4.0D0
         RADIUS2 = 2.5D0**2 + (10.0D0-Y0)**2
         FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS2
      ELSE IF (I.EQ.6) THEN
*
*        lines 20, 21 and 22, belly of letter a, is an ellipse centered
*        in (X0, Y0) with semi-axes 3.5 and 2.75.
*
         X0 = 8.5D0
         Y0 = 5.75D0
         FBND = ((X-X0)/3.5D0)**2 + ((Y-Y0)/2.75D0)**2 - 1.D0
      ELSE IF (I.EQ.7) THEN
*
*        lines 43, 44 and 45, outer curve on bottom of 'g', is an ellipse
*        centered in (X0, Y0) with semi-axes 3.5 and 2.5.
*
         X0 = 17.5D0
         Y0 = 2.5D0
         FBND = ((X-X0)/3.5D0)**2 + ((Y-Y0)/2.5D0)**2 - 1.D0
      ELSE IF (I.EQ.8) THEN
*
*        lines 28, 29 and 30, inner curve on bottom of 'g', is an ellipse
*        centered in (X0, Y0) with semi-axes 2.0 and 1.5.
*
         X0 = 17.5D0
         Y0 = 2.5D0
         FBND = ((X-X0)/2.0D0)**2 + ((Y-Y0)/1.5D0)**2 - 1.D0
      ELSE IF (I.EQ.9) THEN
*
*        line 42, inner curve on lower middle of 'g', is an ellipse
*        centered in (X0, Y0) with semi-axes 1.5 and 0.5.
*
         X0 = 17.5D0
         Y0 = 5.5D0
         FBND = ((X-X0)/1.5D0)**2 + ((Y-Y0)/0.5D0)**2 - 1.D0
      ELSE IF (I.EQ.10) THEN
*
*        line 31, outer curve on lower middle of 'g', is an ellipse
*        centered in (X0, Y0) with semi-axes 2.0 and 1.5.
*
         X0 = 17.5D0
         Y0 = 5.5D0
         FBND = ((X-X0)/3.0D0)**2 + ((Y-Y0)/1.5D0)**2 - 1.D0
      ELSE IF (I.EQ.11) THEN
*
*        line 41, inner curve on upper middle of 'g', is an ellipse
*        centered in (X0, Y0) with semi-axes 1.0 and 1.0.
*
         X0 = 17.0D0
         Y0 = 5.5D0
         FBND = ((X-X0)/1.0D0)**2 + ((Y-Y0)/1.0D0)**2 - 1.D0
      ELSE IF (I.EQ.12) THEN
*
*        line 32, outer curve on upper middle of 'g', is an ellipse
*        centered in (X0, Y0) with semi-axes 1.5 and 1.1573.
*
         X0 = 16.0D0
         Y0 = 5.5D0
         FBND = ((X-X0)/1.5D0)**2 + ((Y-Y0)/1.1573D0)**2 - 1.D0
      ELSE IF (I.EQ.13) THEN
*
*        lines 33, 33, 34, 39 and 40, upper portion of 'g', is an ellipse
*        centered in (X0, Y0) with semi-axes 3.0 and 2.75.
         X0 = 17.0D0
         Y0 = 9.25D0
         FBND = ((X-X0)/3.0D0)**2 + ((Y-Y0)/2.75D0)**2 - 1.D0
      END IF
*
      RETURN
      END