NAG Library Manual, Mark 28.5
```    Program e01eafe

!     E01EAF Example Program Text
!     Mark 28.5 Release. NAG Copyright 2022.

!     .. Use Statements ..
Use nag_library, Only: e01eaf, e01ebf, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: nin = 5, nout = 6
Logical, Parameter               :: pr_tr = .False.
!     .. Local Scalars ..
Integer                          :: i, ifail, j, m, n, nb, nbn, nt
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: f(:), pf(:), px(:), py(:), x(:),     &
y(:)
Integer, Allocatable             :: bnd(:), bnodes(:), tri(:,:),         &
triang(:)
!     .. Executable Statements ..

Write (nout,*) 'E01EAF Example Program Results'

!     Skip heading in data file

Allocate (x(n),y(n),f(n),triang(7*n),tri(3,2*n-5),bnd(2*n-5))

!     Triangulate data
ifail = 0
Call e01eaf(n,x,y,triang,ifail)

!     Convert array triang to an integer array tri(3,2*n-5)
!     such that tri(1:3,k) holds the node indices for the k-th triangle.
!     nt = number of triangles. nb = number of triangles with boundary edge.
!     If bnd(i) > 0 Then tri(:,i) is a triangle with boundary edge;
!       bnd(i) = 4     means that tri(:,i) has two boundary edges,
!       bnd(i) = 1,2,3 means that tri(bnd(i),i) is internal vertex.

Call triang2list(n,triang,nt,tri,nb,bnd)

Write (nout,*)
Write (nout,*) ' Triangle Information'
Write (nout,99999) '  Number of triangles                     = ', nt
Write (nout,99999) '  Number of triangles with boundary edges = ', nb
99999 Format (1X,A,I3)

!     Also find the boundary nodes in anti-clockwise order using triang
!     That is, find the convex hull.
Allocate (bnodes(n))
Call convex_hull(n,triang,bnodes,nbn)

!     Print boundary nodes
Write (nout,*)
Write (nout,*) ' Boundary Node Information'
Write (nout,99999) '  Number of boundary nodes                = ', nbn
Write (nout,*)
Write (nout,*) '  node     Coordinates'
Do i = 1, nbn
j = bnodes(i)
Write (nout,99998) j, x(j), y(j)
End Do
99998 Format (1X,I5,3X,'(',F7.4,', ',F7.4,')')

Allocate (px(m),py(m),pf(m))

!     Interpolate data
ifail = 0
Call e01ebf(m,n,x,y,f,triang,px,py,pf,ifail)

!     Display results
Write (nout,*)
Write (nout,*) ' Interpolation Results'
Write (nout,99997) 'px', 'py', 'Interpolated Value'
Write (nout,99996)(px(i),py(i),pf(i),i=1,m)

If (pr_tr) Then
Call print_triang
End If

99997 Format (2X,A4,4X,A4,4X,A19)
99996 Format (1X,F7.4,1X,F7.4,8X,F7.4)
Contains
Subroutine print_triang

!       .. Implicit None Statement ..
Implicit None
!       .. Local Scalars ..
Integer                        :: i_k, j, j_k, k
!       .. Executable Statements ..

!       Print a sequence of unique line segments for plotting triangulation
Write (nout,*)
Write (nout,*) '  Triangulation as a set of line segments'
Write (nout,*)
j_k = 0
Do k = 1, n
i_k = j_k + 1
j_k = triang(6*n+k)
Do j = i_k, j_k
If (triang(j)>k) Then
Write (nout,99999) x(k), y(k)
Write (nout,99999) x(triang(j)), y(triang(j))
Write (nout,*)
End If
End Do
End Do
Return
99999   Format (1X,F7.4,1X,F7.4)
End Subroutine print_triang
Subroutine triang2list(n,triang,nt,tri,nb,bnd)

!       .. Implicit None Statement ..
Implicit None
!       .. Scalar Arguments ..
Integer, Intent (In)           :: n
Integer, Intent (Out)          :: nb, nt
!       .. Array Arguments ..
Integer, Intent (Out)          :: bnd(2*n-5), tri(3,2*n-5)
Integer, Intent (In)           :: triang(7*n)
!       .. Local Scalars ..
Integer                        :: i, i_k, i_ts, j, jj, j_k, k, m
!       .. Local Arrays ..
Integer                        :: b(3), t(3)
!       .. Executable Statements ..

!       m
!          Initial setting:     0
!          During calculation:  current triangle number
!          On final exit:       total number of unique triangles
m = 0
j_k = 0
nb = 0
Do k = 1, n
!         Which nodes are connected to node k?

!         Node k is connected to nodes i_k to j_k
i_k = j_k + 1
j_k = triang(6*n+k)

!         Let n_k (= j_k - i_k + 1) be the number of nodes

!         first connection setup, first two vertices of first triangle:
!                                 node k and node i_k
t(1) = k
t(2) = triang(i_k)

!         For the remaining connected nodes, we need to know whether node k
!         is internal or on the boundary.
!         An internal node has n_k associated triangles;
!         A boundary node has n_k-1 associated triangles.
If (triang(j_k)>0) Then
!           Node k is internal
!           Connected node order is anti-clockwise;
!           so for the last triangle connected to node k has nodes
!           k, j_k and i_k.
!           We need to loop round final connected point to first point
jj = j_k + 1
Else
!           Node k is a boundary points;
!           There are n_k-1 triangles and no need to look at zero marker.
jj = j_k - 1
End If

!         Now loop over the remaining connecting nodes
Do j = i_k + 1, jj
If (j==j_k+1) Then
!             final triangle; last node is i_k
t(3) = triang(i_k)
Else
t(3) = triang(j)
End If
!           Each triangle will be visited a number of times, but since we
!           are going through in ascending node number k, a new triangle
!           is identified by having connecting node numbers greater than k.
If (t(2)>k .And. t(3)>k) Then
!             new triangle here
m = m + 1
tri(1:3,m) = t(1:3)

!             First assume that no edge of triangle is on boundary
bnd(m) = 0

!             If two if the following terms are zero then
!                the triangle has an edge on the boundary
b(1:3) = triang(6*n+t(1:3))
b(1:3) = triang(b(1:3))
i_ts = b(1) + b(2) + b(3)
If (i_ts==0) Then
!               triangle lies on corner
bnd(m) = 4
nb = nb + 1
Else
Do i = 1, 3
If (b(i)==i_ts) Then
!                   Triangle edge is on boundary
!                   Set bnd(m) to index of vertex that is internal
bnd(m) = i
nb = nb + 1
End If
End Do
End If
End If
!           shuffle down (round, anti-clockwise) for next connected node
!           The last point in current triangle connected to node k becomes
!           the second point in next triangle.
t(2) = t(3)
End Do
End Do
nt = m
Return
End Subroutine triang2list
Subroutine convex_hull(n,triang,bnodes,nbn)

!       .. Implicit None Statement ..
Implicit None
!       .. Scalar Arguments ..
Integer, Intent (In)           :: n
Integer, Intent (Out)          :: nbn
!       .. Array Arguments ..
Integer, Intent (Out)          :: bnodes(n)
Integer, Intent (In)           :: triang(7*n)
!       .. Local Scalars ..
Integer                        :: i_k, j, jj, j_k, k, k_j, m
Logical                        :: found
!       .. Executable Statements ..

!       Algorithm:
!          1. Find first node in sequence that is boundary node: m = 1
!          2. Get list of nodes connecting to current node
!          3. Find first connecting node that is boundary node and
!             make this the current node: m = m + 1
!          4. Repeat steps 2. and 3. until current mode is first node

!       1. Find first boundary node
k = 1
j_k = triang(6*n+k)
Do While (triang(j_k)>0 .And. k<n)
!         Node k is internal, move to next.
k = k + 1
j_k = triang(6*n+k)
End Do

!       Either node k is boundary or something has gone wrong
If (k==n) Then
bnodes(1:nt) = -1
Return
End If

!       We have first boundary node k
m = 1
bnodes(m) = k

!       Loop until we come back round to node k
found = .True.
Do While (found)

!         2. Which nodes are connected to current node?

If (k==1) Then
i_k = 1
Else
i_k = triang(6*n+k-1) + 1
End If
j_k = triang(6*n+k) - 1
!         Node k is connected to nodes triang(i_k) to triang(j_k)

!         3. Find next connecting node that is boundary node

k = 0
Do jj = i_k, j_k
j = triang(jj)
k_j = triang(6*n+j)
If (triang(k_j)==0) Then
k = j
Exit
End If
End Do

!         Check that boundary node found
If (k==0) Then
!           could not find bnode to connect to mode k bnodes(m)
!           set flag in bnodes(m+1) and exit
bnodes(m+1) = -1
found = .False.
Else If (k==bnodes(1)) Then
!           We have gone right round the boundary now, so exit
found = .False.
Else
!           Update new boundary node and cycle
m = m + 1
bnodes(m) = k
End If
End Do
nbn = m
Return
End Subroutine convex_hull

End Program e01eafe
```