= det(N0,N1,N2) >= 0.
!
left = x0 * ( y1 * z2 - y2 * z1 ) &
- y0 * ( x1 * z2 - x2 * z1 ) &
+ z0 * ( x1 * y2 - x2 * y1 ) >= 0.0D+00
return
end
function lstptr ( lpl, nb, list, lptr )
!*****************************************************************************80
!
!! LSTPTR returns the index of NB in the adjacency list.
!
! Discussion:
!
! This function returns the index (LIST pointer) of NB in
! the adjacency list for N0, where LPL = LEND(N0).
!
! This function is identical to the similarly named function in TRIPACK.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer LPL, is LEND(N0).
!
! Input, integer NB, index of the node whose pointer is to
! be returned. NB must be connected to N0.
!
! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), the data
! structure defining the triangulation, created by TRMESH.
!
! Output, integer LSTPTR, pointer such that LIST(LSTPTR) = NB or
! LIST(LSTPTR) = -NB, unless NB is not a neighbor of N0, in which
! case LSTPTR = LPL.
!
! Local parameters:
!
! LP = LIST pointer
! ND = Nodal index
!
implicit none
integer list(*)
integer lp
integer lpl
integer lptr(*)
integer lstptr
integer nb
integer nd
lp = lptr(lpl)
do
nd = list(lp)
if ( nd == nb ) then
exit
end if
lp = lptr(lp)
if ( lp == lpl ) then
exit
end if
end do
lstptr = lp
return
end
function nbcnt ( lpl, lptr )
!*****************************************************************************80
!
!! NBCNT returns the number of neighbors of a node.
!
! Discussion:
!
! This function returns the number of neighbors of a node
! N0 in a triangulation created by TRMESH.
!
! The number of neighbors also gives the order of the Voronoi
! polygon containing the point. Thus, a neighbor count of 6
! means the node is contained in a 6-sided Voronoi region.
!
! This function is identical to the similarly named function in TRIPACK.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer LPL = LIST pointer to the last neighbor of N0;
! LPL = LEND(N0).
!
! Input, integer LPTR(6*(N-2)), pointers associated with LIST.
!
! Output, integer NBCNT, the number of neighbors of N0.
!
! Local parameters:
!
! K = Counter for computing the number of neighbors.
!
! LP = LIST pointer
!
implicit none
integer k
integer lp
integer lpl
integer lptr(*)
integer nbcnt
lp = lpl
k = 1
do
lp = lptr(lp)
if ( lp == lpl ) then
exit
end if
k = k + 1
end do
nbcnt = k
return
end
function nearnd ( p, ist, n, x, y, z, list, lptr, lend, al )
!*****************************************************************************80
!
!! NEARND returns the nearest node to a given point.
!
! Discussion:
!
! Given a point P on the surface of the unit sphere and a
! Delaunay triangulation created by TRMESH, this
! function returns the index of the nearest triangulation
! node to P.
!
! The algorithm consists of implicitly adding P to the
! triangulation, finding the nearest neighbor to P, and
! implicitly deleting P from the triangulation. Thus, it
! is based on the fact that, if P is a node in a Delaunay
! triangulation, the nearest node to P is a neighbor of P.
!
! For large values of N, this procedure will be faster than
! the naive approach of computing the distance from P to every node.
!
! Note that the number of candidates for NEARND (neighbors of P)
! is limited to LMAX defined in the PARAMETER statement below.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, real ( kind = rk ) P(3), the Cartesian coordinates of the point P to
! be located relative to the triangulation. It is assumed
! that P(1)**2 + P(2)**2 + P(3)**2 = 1, that is, that the
! point lies on the unit sphere.
!
! Input, integer IST, the index of the node at which the search
! is to begin. The search time depends on the proximity of this
! node to P. If no good candidate is known, any value between
! 1 and N will do.
!
! Input, integer N, the number of nodes in the triangulation.
! N must be at least 3.
!
! Input, real ( kind = rk ) X(N), Y(N), Z(N), the Cartesian coordinates of
! the nodes.
!
! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
! the data structure defining the triangulation, created by TRMESH.
!
! Output, real ( kind = rk ) AL, the arc length between P and node NEARND.
! Because both points are on the unit sphere, this is also
! the angular separation in radians.
!
! Output, integer NEARND, the index of the nearest node to P.
! NEARND will be 0 if N < 3 or the triangulation data structure
! is invalid.
!
! Local parameters:
!
! B1,B2,B3 = Unnormalized barycentric coordinates returned by TRFIND
! DS1 = (Negative cosine of the) distance from P to N1
! DSR = (Negative cosine of the) distance from P to NR
! DX1,..DZ3 = Components of vectors used by the swap test
! I1,I2,I3 = Nodal indexes of a triangle containing P, or
! the rightmost (I1) and leftmost (I2) visible
! boundary nodes as viewed from P
! L = Length of LISTP/LPTRP and number of neighbors of P
! LMAX = Maximum value of L
! LISTP = Indexes of the neighbors of P
! LPTRP = Array of pointers in 1-1 correspondence with LISTP elements
! LP = LIST pointer to a neighbor of N1 and LISTP pointer
! LP1,LP2 = LISTP indexes (pointers)
! LPL = Pointer to the last neighbor of N1
! N1 = Index of a node visible from P
! N2 = Index of an endpoint of an arc opposite P
! N3 = Index of the node opposite N1->N2
! NN = Local copy of N
! NR = Index of a candidate for the nearest node to P
! NST = Index of the node at which TRFIND begins the search
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer, parameter :: lmax = 25
integer n
real ( kind = rk ) al
real ( kind = rk ) b1
real ( kind = rk ) b2
real ( kind = rk ) b3
real ( kind = rk ) ds1
real ( kind = rk ) dsr
real ( kind = rk ) dx1
real ( kind = rk ) dx2
real ( kind = rk ) dx3
real ( kind = rk ) dy1
real ( kind = rk ) dy2
real ( kind = rk ) dy3
real ( kind = rk ) dz1
real ( kind = rk ) dz2
real ( kind = rk ) dz3
integer i1
integer i2
integer i3
integer ist
integer l
integer lend(n)
integer list(6*(n-2))
integer listp(lmax)
integer lp
integer lp1
integer lp2
integer lpl
integer lptr(6*(n-2))
integer lptrp(lmax)
integer lstptr
integer nearnd
integer n1
integer n2
integer n3
integer nn
integer nr
integer nst
real ( kind = rk ) p(3)
real ( kind = rk ) x(n)
real ( kind = rk ) y(n)
real ( kind = rk ) z(n)
nearnd = 0
al = 0.0D+00
!
! Store local parameters and test for N invalid.
!
nn = n
if ( nn < 3 ) then
return
end if
nst = ist
if ( nst < 1 .or. nn < nst ) then
nst = 1
end if
!
! Find a triangle (I1,I2,I3) containing P, or the rightmost
! (I1) and leftmost (I2) visible boundary nodes as viewed from P.
!
call trfind ( nst, p, n, x, y, z, list, lptr, lend, b1, b2, b3, i1, i2, i3 )
!
! Test for collinear nodes.
!
if ( i1 == 0 ) then
return
end if
!
! Store the linked list of 'neighbors' of P in LISTP and
! LPTRP. I1 is the first neighbor, and 0 is stored as
! the last neighbor if P is not contained in a triangle.
! L is the length of LISTP and LPTRP, and is limited to
! LMAX.
!
if ( i3 /= 0 ) then
listp(1) = i1
lptrp(1) = 2
listp(2) = i2
lptrp(2) = 3
listp(3) = i3
lptrp(3) = 1
l = 3
else
n1 = i1
l = 1
lp1 = 2
listp(l) = n1
lptrp(l) = lp1
!
! Loop on the ordered sequence of visible boundary nodes
! N1 from I1 to I2.
!
do
lpl = lend(n1)
n1 = -list(lpl)
l = lp1
lp1 = l+1
listp(l) = n1
lptrp(l) = lp1
if ( n1 == i2 .or. lmax <= lp1 ) then
exit
end if
end do
l = lp1
listp(l) = 0
lptrp(l) = 1
end if
!
! Initialize variables for a loop on arcs N1-N2 opposite P
! in which new 'neighbors' are 'swapped' in. N1 follows
! N2 as a neighbor of P, and LP1 and LP2 are the LISTP
! indexes of N1 and N2.
!
lp2 = 1
n2 = i1
lp1 = lptrp(1)
n1 = listp(lp1)
!
! Begin loop: find the node N3 opposite N1->N2.
!
do
lp = lstptr ( lend(n1), n2, list, lptr )
if ( 0 <= list(lp) ) then
lp = lptr(lp)
n3 = abs ( list(lp) )
!
! Swap test: Exit the loop if L = LMAX.
!
if ( l == lmax ) then
exit
end if
dx1 = x(n1) - p(1)
dy1 = y(n1) - p(2)
dz1 = z(n1) - p(3)
dx2 = x(n2) - p(1)
dy2 = y(n2) - p(2)
dz2 = z(n2) - p(3)
dx3 = x(n3) - p(1)
dy3 = y(n3) - p(2)
dz3 = z(n3) - p(3)
!
! Swap: Insert N3 following N2 in the adjacency list for P.
! The two new arcs opposite P must be tested.
!
if ( dx3 * ( dy2 * dz1 - dy1 * dz2 ) - &
dy3 * ( dx2 * dz1 - dx1 * dz2 ) + &
dz3 * ( dx2 * dy1 - dx1 * dy2 ) > 0.0D+00 ) then
l = l+1
lptrp(lp2) = l
listp(l) = n3
lptrp(l) = lp1
lp1 = l
n1 = n3
cycle
end if
end if
!
! No swap: Advance to the next arc and test for termination
! on N1 = I1 (LP1 = 1) or N1 followed by 0.
!
if ( lp1 == 1 ) then
exit
end if
lp2 = lp1
n2 = n1
lp1 = lptrp(lp1)
n1 = listp(lp1)
if ( n1 == 0 ) then
exit
end if
end do
!
! Set NR and DSR to the index of the nearest node to P and
! an increasing function (negative cosine) of its distance
! from P, respectively.
!
nr = i1
dsr = -( x(nr) * p(1) + y(nr) * p(2) + z(nr) * p(3) )
do lp = 2, l
n1 = listp(lp)
if ( n1 == 0 ) then
cycle
end if
ds1 = -( x(n1) * p(1) + y(n1) * p(2) + z(n1) * p(3) )
if ( ds1 < dsr ) then
nr = n1
dsr = ds1
end if
end do
dsr = -dsr
dsr = min ( dsr, 1.0D+00 )
al = acos ( dsr )
nearnd = nr
return
end
subroutine optim ( x, y, z, na, list, lptr, lend, nit, iwk, ier )
!*****************************************************************************80
!
!! OPTIM optimizes the quadrilateral portion of a triangulation.
!
! Discussion:
!
! Given a set of NA triangulation arcs, this subroutine
! optimizes the portion of the triangulation consisting of
! the quadrilaterals (pairs of adjacent triangles) which
! have the arcs as diagonals by applying the circumcircle
! test and appropriate swaps to the arcs.
!
! An iteration consists of applying the swap test and
! swaps to all NA arcs in the order in which they are
! stored. The iteration is repeated until no swap occurs
! or NIT iterations have been performed. The bound on the
! number of iterations may be necessary to prevent an
! infinite loop caused by cycling (reversing the effect of a
! previous swap) due to floating point inaccuracy when four
! or more nodes are nearly cocircular.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, real ( kind = rk ) X(*), Y(*), Z(*), the nodal coordinates.
!
! Input, integer NA, the number of arcs in the set. NA >= 0.
!
! Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
! the data structure defining the triangulation, created by TRMESH.
! On output, updated to reflect the swaps.
!
! Input/output, integer NIT. On input, the maximum number of
! iterations to be performed. NIT = 4*NA should be sufficient. NIT >= 1.
! On output, the number of iterations performed.
!
! Input/output, integer IWK(2,NA), the nodal indexes of the arc
! endpoints (pairs of endpoints are stored in columns). On output, endpoint
! indexes of the new set of arcs reflecting the swaps.
!
! Output, integer IER, error indicator:
! 0, if no errors were encountered.
! 1, if a swap occurred on the last of MAXIT iterations, where MAXIT is the
! value of NIT on input. The new set of arcs is not necessarily optimal
! in this case.
! 2, if NA < 0 or NIT < 1 on input.
! 3, if IWK(2,I) is not a neighbor of IWK(1,I) for some I in the range 1
! to NA. A swap may have occurred in this case.
! 4, if a zero pointer was returned by subroutine SWAP.
!
! Local parameters:
!
! I = Column index for IWK
! IO1,IO2 = Nodal indexes of the endpoints of an arc in IWK
! ITER = Iteration count
! LP = LIST pointer
! LP21 = Parameter returned by SWAP (not used)
! LPL = Pointer to the last neighbor of IO1
! LPP = Pointer to the node preceding IO2 as a neighbor of IO1
! MAXIT = Input value of NIT
! N1,N2 = Nodes opposite IO1->IO2 and IO2->IO1, respectively
! NNA = Local copy of NA
! SWP = Flag set to TRUE iff a swap occurs in the optimization loop
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer na
integer i
integer ier
integer io1
integer io2
integer iter
integer iwk(2,na)
integer lend(*)
integer list(*)
integer lp
integer lp21
integer lpl
integer lpp
integer lptr(*)
integer maxit
integer n1
integer n2
integer nit
integer nna
logical swp
logical swptst
real ( kind = rk ) x(*)
real ( kind = rk ) y(*)
real ( kind = rk ) z(*)
nna = na
maxit = nit
if ( nna < 0 .or. maxit < 1 ) then
nit = 0
ier = 2
return
end if
!
! Initialize iteration count ITER and test for NA = 0.
!
iter = 0
if ( nna == 0 ) then
nit = 0
ier = 0
return
end if
!
! Top of loop.
! SWP = TRUE iff a swap occurred in the current iteration.
!
do
if ( maxit <= iter ) then
nit = iter
ier = 1
return
end if
iter = iter + 1
swp = .false.
!
! Inner loop on arcs IO1-IO2.
!
do i = 1, nna
io1 = iwk(1,i)
io2 = iwk(2,i)
!
! Set N1 and N2 to the nodes opposite IO1->IO2 and
! IO2->IO1, respectively. Determine the following:
!
! LPL = pointer to the last neighbor of IO1,
! LP = pointer to IO2 as a neighbor of IO1, and
! LPP = pointer to the node N2 preceding IO2.
!
lpl = lend(io1)
lpp = lpl
lp = lptr(lpp)
do
if ( list(lp) == io2 ) then
go to 3
end if
lpp = lp
lp = lptr(lpp)
if ( lp == lpl ) then
exit
end if
end do
!
! IO2 should be the last neighbor of IO1. Test for no
! arc and bypass the swap test if IO1 is a boundary node.
!
if ( abs ( list(lp) ) /= io2 ) then
nit = iter
ier = 3
return
end if
if ( list(lp) < 0 ) then
go to 4
end if
!
! Store N1 and N2, or bypass the swap test if IO1 is a
! boundary node and IO2 is its first neighbor.
!
3 continue
n2 = list(lpp)
!
! Test IO1-IO2 for a swap, and update IWK if necessary.
!
if ( 0 <= n2 ) then
lp = lptr(lp)
n1 = abs ( list(lp) )
if ( swptst ( n1, n2, io1, io2, x, y, z ) ) then
call swap ( n1, n2, io1, io2, list, lptr, lend, lp21 )
if ( lp21 == 0 ) then
nit = iter
ier = 4
return
end if
swp = .true.
iwk(1,i) = n1
iwk(2,i) = n2
end if
end if
4 continue
end do
if ( .not. swp ) then
exit
end if
end do
nit = iter
ier = 0
return
end
subroutine r83vec_normalize ( n, x, y, z )
!*****************************************************************************80
!
!! R83VEC_NORMALIZE normalizes each R83 in an R83VEC to have unit norm.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 25 June 2002
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of vectors.
!
! Input/output, real ( kind = rk ) X(N), Y(N), Z(N), the components of
! the vectors.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer n
integer i
real ( kind = rk ) norm
real ( kind = rk ) x(n)
real ( kind = rk ) y(n)
real ( kind = rk ) z(n)
do i = 1, n
norm = sqrt ( x(i)**2 + y(i)**2 + z(i)**2 )
if ( norm /= 0.0D+00 ) then
x(i) = x(i) / norm
y(i) = y(i) / norm
z(i) = z(i) / norm
end if
end do
return
end
subroutine scoord ( px, py, pz, plat, plon, pnrm )
!*****************************************************************************80
!
!! SCOORD converts from Cartesian to spherical coordinates.
!
! Discussion:
!
! This subroutine converts a point P from Cartesian (X,Y,Z) coordinates
! to spherical ( LATITUDE, LONGITUDE, RADIUS ) coordinates.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, real ( kind = rk ) PX, PY, PZ, the coordinates of P.
!
! Output, real ( kind = rk ) PLAT, the latitude of P in the range -PI/2
! to PI/2, or 0 if PNRM = 0.
!
! Output, real ( kind = rk ) PLON, the longitude of P in the range -PI to PI,
! or 0 if P lies on the Z-axis.
!
! Output, real ( kind = rk ) PNRM, the magnitude (Euclidean norm) of P.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) plat
real ( kind = rk ) plon
real ( kind = rk ) pnrm
real ( kind = rk ) px
real ( kind = rk ) py
real ( kind = rk ) pz
pnrm = sqrt ( px * px + py * py + pz * pz )
if ( px /= 0.0D+00 .or. py /= 0.0D+00 ) then
plon = atan2 ( py, px )
else
plon = 0.0D+00
end if
if ( pnrm /= 0.0D+00 ) then
plat = asin ( pz / pnrm )
else
plat = 0.0D+00
end if
return
end
function store ( x )
!*****************************************************************************80
!
!! STORE forces its argument to be stored.
!
! Discussion:
!
! This function forces its argument X to be stored in a
! memory location, thus providing a means of determining
! floating point number characteristics (such as the machine
! precision) when it is necessary to avoid computation in
! high precision registers.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, real ( kind = rk ) X, the value to be stored.
!
! Output, real ( kind = rk ) STORE, the value of X after it has been stored
! and possibly truncated or rounded to the single precision word length.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) store
real ( kind = rk ) x
real ( kind = rk ) y
common /stcom/ y
y = x
store = y
return
end
subroutine swap ( in1, in2, io1, io2, list, lptr, lend, lp21 )
!*****************************************************************************80
!
!! SWAP replaces the diagonal arc of a quadrilateral with the other diagonal.
!
! Discussion:
!
! Given a triangulation of a set of points on the unit
! sphere, this subroutine replaces a diagonal arc in a
! strictly convex quadrilateral (defined by a pair of adja-
! cent triangles) with the other diagonal. Equivalently, a
! pair of adjacent triangles is replaced by another pair
! having the same union.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer IN1, IN2, IO1, IO2, nodal indexes of the
! vertices of the quadrilateral. IO1-IO2 is replaced by IN1-IN2.
! (IO1,IO2,IN1) and (IO2,IO1,IN2) must be triangles on input.
!
! Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
! the data structure defining the triangulation, created by TRMESH.
! On output, updated with the swap; triangles (IO1,IO2,IN1) an (IO2,IO1,IN2)
! are replaced by (IN1,IN2,IO2) and (IN2,IN1,IO1) unless LP21 = 0.
!
! Output, integer LP21, index of IN1 as a neighbor of IN2 after
! the swap is performed unless IN1 and IN2 are adjacent on input, in which
! case LP21 = 0.
!
! Local parameters:
!
! LP, LPH, LPSAV = LIST pointers
!
implicit none
integer in1
integer in2
integer io1
integer io2
integer lend(*)
integer list(*)
integer lp
integer lp21
integer lph
integer lpsav
integer lptr(*)
integer lstptr
!
! Test for IN1 and IN2 adjacent.
!
lp = lstptr ( lend(in1), in2, list, lptr )
if ( abs ( list(lp) ) == in2 ) then
lp21 = 0
return
end if
!
! Delete IO2 as a neighbor of IO1.
!
lp = lstptr ( lend(io1), in2, list, lptr )
lph = lptr(lp)
lptr(lp) = lptr(lph)
!
! If IO2 is the last neighbor of IO1, make IN2 the last neighbor.
!
if ( lend(io1) == lph ) then
lend(io1) = lp
end if
!
! Insert IN2 as a neighbor of IN1 following IO1 using the hole created above.
!
lp = lstptr ( lend(in1), io1, list, lptr )
lpsav = lptr(lp)
lptr(lp) = lph
list(lph) = in2
lptr(lph) = lpsav
!
! Delete IO1 as a neighbor of IO2.
!
lp = lstptr ( lend(io2), in1, list, lptr )
lph = lptr(lp)
lptr(lp) = lptr(lph)
!
! If IO1 is the last neighbor of IO2, make IN1 the last neighbor.
!
if ( lend(io2) == lph ) then
lend(io2) = lp
end if
!
! Insert IN1 as a neighbor of IN2 following IO2.
!
lp = lstptr ( lend(in2), io2, list, lptr )
lpsav = lptr(lp)
lptr(lp) = lph
list(lph) = in1
lptr(lph) = lpsav
lp21 = lph
return
end
function swptst ( n1, n2, n3, n4, x, y, z )
!*****************************************************************************80
!
!! SWPTST decides whether to replace a diagonal arc by the other.
!
! Discussion:
!
! This function decides whether or not to replace a
! diagonal arc in a quadrilateral with the other diagonal.
! The decision will be to swap (SWPTST = TRUE) if and only
! if N4 lies above the plane (in the half-space not contain-
! ing the origin) defined by (N1,N2,N3), or equivalently, if
! the projection of N4 onto this plane is interior to the
! circumcircle of (N1,N2,N3). The decision will be for no
! swap if the quadrilateral is not strictly convex.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer N1, N2, N3, N4, indexes of the four nodes
! defining the quadrilateral with N1 adjacent to N2, and (N1,N2,N3) in
! counterclockwise order. The arc connecting N1 to N2 should be replaced
! by an arc connecting N3 to N4 if SWPTST = TRUE. Refer to subroutine SWAP.
!
! Input, real ( kind = rk ) X(N), Y(N), Z(N), the coordinates of the nodes.
!
! Output, logical SWPTST, TRUE if and only if the arc connecting N1
! and N2 should be swapped for an arc connecting N3 and N4.
!
! Local parameters:
!
! DX1,DY1,DZ1 = Coordinates of N4->N1
! DX2,DY2,DZ2 = Coordinates of N4->N2
! DX3,DY3,DZ3 = Coordinates of N4->N3
! X4,Y4,Z4 = Coordinates of N4
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) dx1
real ( kind = rk ) dx2
real ( kind = rk ) dx3
real ( kind = rk ) dy1
real ( kind = rk ) dy2
real ( kind = rk ) dy3
real ( kind = rk ) dz1
real ( kind = rk ) dz2
real ( kind = rk ) dz3
integer n1
integer n2
integer n3
integer n4
logical swptst
real ( kind = rk ) x(*)
real ( kind = rk ) x4
real ( kind = rk ) y(*)
real ( kind = rk ) y4
real ( kind = rk ) z(*)
real ( kind = rk ) z4
x4 = x(n4)
y4 = y(n4)
z4 = z(n4)
dx1 = x(n1) - x4
dx2 = x(n2) - x4
dx3 = x(n3) - x4
dy1 = y(n1) - y4
dy2 = y(n2) - y4
dy3 = y(n3) - y4
dz1 = z(n1) - z4
dz2 = z(n2) - z4
dz3 = z(n3) - z4
!
! N4 lies above the plane of (N1,N2,N3) iff N3 lies above
! the plane of (N2,N1,N4) iff Det(N3-N4,N2-N4,N1-N4) =
! (N3-N4,N2-N4 X N1-N4) > 0.
!
swptst = dx3 * ( dy2 * dz1 - dy1 * dz2 ) &
- dy3 * ( dx2 * dz1 - dx1 * dz2 ) &
+ dz3 * ( dx2 * dy1 - dx1 * dy2 ) > 0.0D+00
return
end
subroutine timestamp ( )
!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
! Example:
!
! 31 May 2001 9:45:54.872 AM
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 18 May 2013
!
! Author:
!
! John Burkardt
!
implicit none
character ( len = 8 ) ampm
integer d
integer h
integer m
integer mm
character ( len = 9 ), parameter, dimension(12) :: month = (/ &
'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', &
'September', 'October ', 'November ', 'December ' /)
integer n
integer s
integer values(8)
integer y
call date_and_time ( values = values )
y = values(1)
m = values(2)
d = values(3)
h = values(5)
n = values(6)
s = values(7)
mm = values(8)
if ( h < 12 ) then
ampm = 'AM'
else if ( h == 12 ) then
if ( n == 0 .and. s == 0 ) then
ampm = 'Noon'
else
ampm = 'PM'
end if
else
h = h - 12
if ( h < 12 ) then
ampm = 'PM'
else if ( h == 12 ) then
if ( n == 0 .and. s == 0 ) then
ampm = 'Midnight'
else
ampm = 'AM'
end if
end if
end if
write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )
return
end
subroutine trans ( n, rlat, rlon, x, y, z )
!*****************************************************************************80
!
!! TRANS transforms spherical coordinates to Cartesian coordinates.
!
! Discussion:
!
! This subroutine transforms spherical coordinates into
! Cartesian coordinates on the unit sphere for input to
! TRMESH. Storage for X and Y may coincide with
! storage for RLAT and RLON if the latter need not be saved.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer N, the number of nodes (points on the unit
! sphere) whose coordinates are to be transformed.
!
! Input, real ( kind = rk ) RLAT(N), latitudes of the nodes in radians.
!
! Input, real ( kind = rk ) RLON(N), longitudes of the nodes in radians.
!
! Output, real ( kind = rk ) X(N), Y(N), Z(N), the coordinates in the
! range -1 to 1. X(I)**2 + Y(I)**2 + Z(I)**2 = 1 for I = 1 to N.
!
! Local parameters:
!
! COSPHI = cos(PHI)
! I = DO-loop index
! NN = Local copy of N
! PHI = Latitude
! THETA = Longitude
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer n
real ( kind = rk ) cosphi
integer i
integer nn
real ( kind = rk ) phi
real ( kind = rk ) rlat(n)
real ( kind = rk ) rlon(n)
real ( kind = rk ) theta
real ( kind = rk ) x(n)
real ( kind = rk ) y(n)
real ( kind = rk ) z(n)
nn = n
do i = 1, nn
phi = rlat(i)
theta = rlon(i)
cosphi = cos ( phi )
x(i) = cosphi * cos ( theta )
y(i) = cosphi * sin ( theta )
z(i) = sin ( phi )
end do
return
end
subroutine trfind ( nst, p, n, x, y, z, list, lptr, lend, b1, b2, b3, i1, &
i2, i3 )
!*****************************************************************************80
!
!! TRFIND locates a point relative to a triangulation.
!
! Discussion:
!
! This subroutine locates a point P relative to a triangulation
! created by TRMESH. If P is contained in
! a triangle, the three vertex indexes and barycentric
! coordinates are returned. Otherwise, the indexes of the
! visible boundary nodes are returned.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer NST, index of a node at which TRFIND begins
! its search. Search time depends on the proximity of this node to P.
!
! Input, real ( kind = rk ) P(3), the x, y, and z coordinates (in that order)
! of the point P to be located.
!
! Input, integer N, the number of nodes in the triangulation.
! 3 <= N.
!
! Input, real ( kind = rk ) X(N), Y(N), Z(N), the coordinates of the
! triangulation nodes (unit vectors).
!
! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), the
! data structure defining the triangulation, created by TRMESH.
!
! Output, real ( kind = rk ) B1, B2, B3, the unnormalized barycentric
! coordinates of the central projection of P onto the underlying planar
! triangle if P is in the convex hull of the nodes. These parameters
! are not altered if I1 = 0.
!
! Output, integer I1, I2, I3, the counterclockwise-ordered
! vertex indexes of a triangle containing P if P is contained in a triangle.
! If P is not in the convex hull of the nodes, I1 and I2 are the rightmost
! and leftmost (boundary) nodes that are visible from P, and I3 = 0. (If
! all boundary nodes are visible from P, then I1 and I2 coincide.)
! I1 = I2 = I3 = 0 if P and all of the nodes are coplanar (lie on a
! common great circle.
!
! Local parameters:
!
! EPS = Machine precision
! IX,IY,IZ = Integer seeds for JRAND
! LP = LIST pointer
! N0,N1,N2 = Nodes in counterclockwise order defining a
! cone (with vertex N0) containing P, or end-
! points of a boundary edge such that P Right
! N1->N2
! N1S,N2S = Initially-determined values of N1 and N2
! N3,N4 = Nodes opposite N1->N2 and N2->N1, respectively
! NEXT = Candidate for I1 or I2 when P is exterior
! NF,NL = First and last neighbors of N0, or first
! (rightmost) and last (leftmost) nodes
! visible from P when P is exterior to the
! triangulation
! PTN1 = Scalar product
! PTN2 = Scalar product

! Q = (N2 X N1) X N2 or N1 X (N2 X N1) -- used in
! the boundary traversal when P is exterior
! S12 = Scalar product
! TOL = Tolerance (multiple of EPS) defining an upper
! bound on the magnitude of a negative bary-
! centric coordinate (B1 or B2) for P in a
! triangle -- used to avoid an infinite number
! of restarts with 0 <= B3 < EPS and B1 < 0 or
! B2 < 0 but small in magnitude
! XP,YP,ZP = Local variables containing P(1), P(2), and P(3)
! X0,Y0,Z0 = Dummy arguments for DET
! X1,Y1,Z1 = Dummy arguments for DET
! X2,Y2,Z2 = Dummy arguments for DET
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer n
real ( kind = rk ) b1
real ( kind = rk ) b2
real ( kind = rk ) b3
real ( kind = rk ) det
real ( kind = rk ) eps
integer i1
integer i2
integer i3
integer, save :: ix = 1
integer, save :: iy = 2
integer, save :: iz = 3
integer jrand
integer lend(n)
integer list(6*(n-2))
integer lp
integer lptr(6*(n-2))
integer lstptr
integer n0
integer n1
integer n1s
integer n2
integer n2s
integer n3
integer n4
integer next
integer nf
integer nl
integer nst
real ( kind = rk ) p(3)
real ( kind = rk ) ptn1
real ( kind = rk ) ptn2
real ( kind = rk ) q(3)
real ( kind = rk ) s12
real ( kind = rk ) store
real ( kind = rk ) tol
real ( kind = rk ) x(n)
real ( kind = rk ) x0
real ( kind = rk ) x1
real ( kind = rk ) x2
real ( kind = rk ) xp
real ( kind = rk ) y(n)
real ( kind = rk ) y0
real ( kind = rk ) y1
real ( kind = rk ) y2
real ( kind = rk ) yp
real ( kind = rk ) z(n)
real ( kind = rk ) z0
real ( kind = rk ) z1
real ( kind = rk ) z2
real ( kind = rk ) zp
!
! Statement function:
!
! DET(X1,...,Z0) >= 0 if and only if (X0,Y0,Z0) is in the
! (closed) left hemisphere defined by the plane containing (0,0,0),
! (X1,Y1,Z1), and (X2,Y2,Z2), where left is defined relative to an
! observer at (X1,Y1,Z1) facing (X2,Y2,Z2).
!
det (x1,y1,z1,x2,y2,z2,x0,y0,z0) = x0*(y1*z2-y2*z1) &
- y0*(x1*z2-x2*z1) + z0*(x1*y2-x2*y1)
!
! Initialize variables.
!
xp = p(1)
yp = p(2)
zp = p(3)
n0 = nst
if ( n0 < 1 .or. n < n0 ) then
n0 = jrand ( n, ix, iy, iz )
end if
!
! Compute the relative machine precision EPS and TOL.
!
eps = epsilon ( eps )
tol = 100.0D+00 * eps
!
! Set NF and NL to the first and last neighbors of N0, and initialize N1 = NF.
!
2 continue
lp = lend(n0)
nl = list(lp)
lp = lptr(lp)
nf = list(lp)
n1 = nf
!
! Find a pair of adjacent neighbors N1,N2 of N0 that define
! a wedge containing P: P LEFT N0->N1 and P RIGHT N0->N2.
!
if ( 0 < nl ) then
!
! N0 is an interior node. Find N1.
!
3 continue
if ( det ( x(n0),y(n0),z(n0),x(n1),y(n1),z(n1),xp,yp,zp ) < 0.0D+00 ) then
lp = lptr(lp)
n1 = list(lp)
if ( n1 == nl ) then
go to 6
end if
go to 3
end if
else
!
! N0 is a boundary node. Test for P exterior.
!
nl = -nl
!
! Is P to the right of the boundary edge N0->NF?
!
if ( det(x(n0),y(n0),z(n0),x(nf),y(nf),z(nf), xp,yp,zp) < 0.0D+00 ) then
n1 = n0
n2 = nf
go to 9
end if
!
! Is P to the right of the boundary edge NL->N0?
!
if ( det(x(nl),y(nl),z(nl),x(n0),y(n0),z(n0),xp,yp,zp) < 0.0D+00 ) then
n1 = nl
n2 = n0
go to 9
end if
end if
!
! P is to the left of arcs N0->N1 and NL->N0. Set N2 to the
! next neighbor of N0 (following N1).
!
4 continue
lp = lptr(lp)
n2 = abs ( list(lp) )
if ( det(x(n0),y(n0),z(n0),x(n2),y(n2),z(n2),xp,yp,zp) < 0.0D+00 ) then
go to 7
end if
n1 = n2
if ( n1 /= nl ) then
go to 4
end if
if ( det ( x(n0), y(n0), z(n0), x(nf), y(nf), z(nf), xp, yp, zp ) &
< 0.0D+00 ) then
go to 6
end if
!
! P is left of or on arcs N0->NB for all neighbors NB
! of N0. Test for P = +/-N0.
!
if ( store ( abs ( x(n0 ) * xp + y(n0) * yp + z(n0) * zp) ) &
< 1.0D+00 - 4.0D+00 * eps ) then
!
! All points are collinear iff P Left NB->N0 for all
! neighbors NB of N0. Search the neighbors of N0.
! Note: N1 = NL and LP points to NL.
!
do
if ( det(x(n1),y(n1),z(n1),x(n0),y(n0),z(n0),xp,yp,zp) < 0.0D+00 ) then
exit
end if
lp = lptr(lp)
n1 = abs ( list(lp) )
if ( n1 == nl ) then
i1 = 0
i2 = 0
i3 = 0
return
end if
end do
end if
!
! P is to the right of N1->N0, or P = +/-N0. Set N0 to N1 and start over.
!
n0 = n1
go to 2
!
! P is between arcs N0->N1 and N0->NF.
!
6 continue
n2 = nf
!
! P is contained in a wedge defined by geodesics N0-N1 and
! N0-N2, where N1 is adjacent to N2. Save N1 and N2 to
! test for cycling.
!
7 continue
n3 = n0
n1s = n1
n2s = n2
!
! Top of edge-hopping loop:
!
8 continue
b3 = det ( x(n1),y(n1),z(n1),x(n2),y(n2),z(n2),xp,yp,zp )
if ( b3 < 0.0D+00 ) then
!
! Set N4 to the first neighbor of N2 following N1 (the
! node opposite N2->N1) unless N1->N2 is a boundary arc.
!
lp = lstptr ( lend(n2), n1, list, lptr )
if ( list(lp) < 0 ) then
go to 9
end if
lp = lptr(lp)
n4 = abs ( list(lp) )
!
! Define a new arc N1->N2 which intersects the geodesic N0-P.
!
if ( det ( x(n0),y(n0),z(n0),x(n4),y(n4),z(n4),xp,yp,zp ) < 0.0D+00 ) then
n3 = n2
n2 = n4
n1s = n1
if ( n2 /= n2s .and. n2 /= n0 ) then
go to 8
end if
else
n3 = n1
n1 = n4
n2s = n2
if ( n1 /= n1s .and. n1 /= n0 ) then
go to 8
end if
end if
!
! The starting node N0 or edge N1-N2 was encountered
! again, implying a cycle (infinite loop). Restart
! with N0 randomly selected.
!
n0 = jrand ( n, ix, iy, iz )
go to 2
end if
!
! P is in (N1,N2,N3) unless N0, N1, N2, and P are collinear
! or P is close to -N0.
!
if ( b3 >= eps ) then
!
! B3 /= 0.
!
b1 = det(x(n2),y(n2),z(n2),x(n3),y(n3),z(n3),xp,yp,zp)
b2 = det(x(n3),y(n3),z(n3),x(n1),y(n1),z(n1),xp,yp,zp)
!
! Restart with N0 randomly selected.
!
if ( b1 < -tol .or. b2 < -tol ) then
n0 = jrand ( n, ix, iy, iz )
go to 2
end if
else
!
! B3 = 0 and thus P lies on N1->N2. Compute
! B1 = Det(P,N2 X N1,N2) and B2 = Det(P,N1,N2 X N1).
!
b3 = 0.0D+00
s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
ptn1 = xp * x(n1) + yp * y(n1) + zp * z(n1)
ptn2 = xp * x(n2) + yp * y(n2) + zp * z(n2)
b1 = ptn1 - s12 * ptn2
b2 = ptn2 - s12 * ptn1
!
! Restart with N0 randomly selected.
!
if ( b1 < -tol .or. b2 < -tol ) then
n0 = jrand ( n, ix, iy, iz )
go to 2
end if
end if
!
! P is in (N1,N2,N3).
!
i1 = n1
i2 = n2
i3 = n3
b1 = max ( b1, 0.0D+00 )
b2 = max ( b2, 0.0D+00 )
return
!
! P Right N1->N2, where N1->N2 is a boundary edge.
! Save N1 and N2, and set NL = 0 to indicate that
! NL has not yet been found.
!
9 continue
n1s = n1
n2s = n2
nl = 0
!
! Counterclockwise Boundary Traversal:
!
10 continue
lp = lend(n2)
lp = lptr(lp)
next = list(lp)
if ( det(x(n2),y(n2),z(n2),x(next),y(next),z(next),xp,yp,zp) >= 0.0D+00 ) then
!
! N2 is the rightmost visible node if P Forward N2->N1
! or NEXT Forward N2->N1. Set Q to (N2 X N1) X N2.
!
s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
q(1) = x(n1) - s12 * x(n2)
q(2) = y(n1) - s12 * y(n2)
q(3) = z(n1) - s12 * z(n2)
if ( xp * q(1) + yp * q(2) + zp * q(3) >= 0.0D+00 ) then
go to 11
end if
if ( x(next) * q(1) + y(next) * q(2) + z(next) * q(3) >= 0.0D+00 ) then
go to 11
end if
!
! N1, N2, NEXT, and P are nearly collinear, and N2 is
! the leftmost visible node.
!
nl = n2
end if
!
! Bottom of counterclockwise loop:
!
n1 = n2
n2 = next
if ( n2 /= n1s ) then
go to 10
end if
!
! All boundary nodes are visible from P.
!
i1 = n1s
i2 = n1s
i3 = 0
return
!
! N2 is the rightmost visible node.
!
11 continue
nf = n2
if ( nl == 0 ) then
!
! Restore initial values of N1 and N2, and begin the search
! for the leftmost visible node.
!
n2 = n2s
n1 = n1s
!
! Clockwise Boundary Traversal:
!
12 continue
lp = lend(n1)
next = -list(lp)
if ( 0.0D+00 <= &
det ( x(next), y(next), z(next), x(n1), y(n1), z(n1), xp, yp, zp ) ) then
!
! N1 is the leftmost visible node if P or NEXT is
! forward of N1->N2. Compute Q = N1 X (N2 X N1).
!
s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
q(1) = x(n2) - s12 * x(n1)
q(2) = y(n2) - s12 * y(n1)
q(3) = z(n2) - s12 * z(n1)
if ( xp * q(1) + yp * q(2) + zp * q(3) >= 0.0D+00 ) then
go to 13
end if
if ( x(next) * q(1) + y(next) * q(2) + z(next) * q(3) >= 0.0D+00 ) then
go to 13
end if
!
! P, NEXT, N1, and N2 are nearly collinear and N1 is the rightmost
! visible node.
!
nf = n1
end if
!
! Bottom of clockwise loop:
!
n2 = n1
n1 = next
if ( n1 /= n1s ) then
go to 12
end if
!
! All boundary nodes are visible from P.
!
i1 = n1
i2 = n1
i3 = 0
return
!
! N1 is the leftmost visible node.
!
13 continue
nl = n1
end if
!
! NF and NL have been found.
!
i1 = nf
i2 = nl
i3 = 0
return
end
subroutine trlist ( n, list, lptr, lend, nrow, nt, ltri, ier )
!*****************************************************************************80
!
!! TRLIST converts a triangulation data structure to a triangle list.
!
! Discussion:
!
! This subroutine converts a triangulation data structure
! from the linked list created by TRMESH to a triangle list.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer N, the number of nodes in the triangulation.
! 3 <= N.
!
! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), linked
! list data structure defining the triangulation. Refer to TRMESH.
!
! Input, integer NROW, the number of rows (entries per triangle)
! reserved for the triangle list LTRI. The value must be 6 if only the
! vertex indexes and neighboring triangle indexes are to be stored, or 9
! if arc indexes are also to be assigned and stored. Refer to LTRI.
!
! Output, integer NT, the number of triangles in the
! triangulation unless IER /=0, in which case NT = 0. NT = 2N-NB-2 if
! NB >= 3 or 2N-4 if NB = 0, where NB is the number of boundary nodes.
!
! Output, integer LTRI(NROW,*). The second dimension of LTRI
! must be at least NT, where NT will be at most 2*N-4. The J-th column
! contains the vertex nodal indexes (first three rows), neighboring triangle
! indexes (second three rows), and, if NROW = 9, arc indexes (last three
! rows) associated with triangle J for J = 1,...,NT. The vertices are
! ordered counterclockwise with the first vertex taken to be the one
! with smallest index. Thus, LTRI(2,J) and LTRI(3,J) are larger than
! LTRI(1,J) and index adjacent neighbors of node LTRI(1,J). For
! I = 1,2,3, LTRI(I+3,J) and LTRI(I+6,J) index the triangle and arc,
! respectively, which are opposite (not shared by) node LTRI(I,J), with
! LTRI(I+3,J) = 0 if LTRI(I+6,J) indexes a boundary arc. Vertex indexes
! range from 1 to N, triangle indexes from 0 to NT, and, if included,
! arc indexes from 1 to NA, where NA = 3N-NB-3 if NB >= 3 or 3N-6 if
! NB = 0. The triangles are ordered on first (smallest) vertex indexes.
!
! Output, integer IER, error indicator.
! 0, if no errors were encountered.
! 1, if N or NROW is outside its valid range on input.
! 2, if the triangulation data structure (LIST,LPTR,LEND) is invalid.
! Note, however, that these arrays are not completely tested for validity.
!
! Local parameters:
!
! ARCS = Logical variable with value TRUE iff are
! indexes are to be stored
! I,J = LTRI row indexes (1 to 3) associated with
! triangles KT and KN, respectively
! I1,I2,I3 = Nodal indexes of triangle KN
! ISV = Variable used to permute indexes I1,I2,I3
! KA = Arc index and number of currently stored arcs
! KN = Index of the triangle that shares arc I1-I2 with KT
! KT = Triangle index and number of currently stored triangles
! LP = LIST pointer
! LP2 = Pointer to N2 as a neighbor of N1
! LPL = Pointer to the last neighbor of I1
! LPLN1 = Pointer to the last neighbor of N1
! N1,N2,N3 = Nodal indexes of triangle KT
! NM2 = N-2
!
implicit none
integer n
integer nrow
logical arcs
integer i
integer i1
integer i2
integer i3
integer ier
integer isv
integer j
integer ka
integer kn
integer kt
integer lend(n)
integer list(6*(n-2))
integer lp
integer lp2
integer lpl
integer lpln1
integer lptr(6*(n-2))
integer ltri(nrow,*)
integer n1
integer n2
integer n3
integer nm2
integer nt
!
! Test for invalid input parameters.
!
if ( n < 3 .or. ( nrow /= 6 .and. nrow /= 9 ) ) then
nt = 0
ier = 1
return
end if
!
! Initialize parameters for loop on triangles KT = (N1,N2,
! N3), where N1 < N2 and N1 < N3.
!
! ARCS = TRUE iff arc indexes are to be stored.
! KA,KT = Numbers of currently stored arcs and triangles.
! NM2 = Upper bound on candidates for N1.
!
arcs = nrow == 9
ka = 0
kt = 0
nm2 = n-2
!
! Loop on nodes N1.
!
do n1 = 1, nm2
!
! Loop on pairs of adjacent neighbors (N2,N3). LPLN1 points
! to the last neighbor of N1, and LP2 points to N2.
!
lpln1 = lend(n1)
lp2 = lpln1
1 continue
lp2 = lptr(lp2)
n2 = list(lp2)
lp = lptr(lp2)
n3 = abs ( list(lp) )
if ( n2 < n1 .or. n3 < n1 ) then
go to 8
end if
!
! Add a new triangle KT = (N1,N2,N3).
!
kt = kt + 1
ltri(1,kt) = n1
ltri(2,kt) = n2
ltri(3,kt) = n3
!
! Loop on triangle sides (I2,I1) with neighboring triangles
! KN = (I1,I2,I3).
!
do i = 1, 3
if ( i == 1 ) then
i1 = n3
i2 = n2
else if ( i == 2 ) then
i1 = n1
i2 = n3
else
i1 = n2
i2 = n1
end if
!
! Set I3 to the neighbor of I1 that follows I2 unless
! I2->I1 is a boundary arc.
!
lpl = lend(i1)
lp = lptr(lpl)
do
if ( list(lp) == i2 ) then
go to 3
end if
lp = lptr(lp)
if ( lp == lpl ) then
exit
end if
end do
!
! Invalid triangulation data structure: I1 is a neighbor of
! I2, but I2 is not a neighbor of I1.
!
if ( abs ( list(lp) ) /= i2 ) then
nt = 0
ier = 2
return
end if
!
! I2 is the last neighbor of I1. Bypass the search for a neighboring
! triangle if I2->I1 is a boundary arc.
!
kn = 0
if ( list(lp) < 0 ) then
go to 6
end if
!
! I2->I1 is not a boundary arc, and LP points to I2 as
! a neighbor of I1.
!
3 continue
lp = lptr(lp)
i3 = abs ( list(lp) )
!
! Find J such that LTRI(J,KN) = I3 (not used if KT < KN),
! and permute the vertex indexes of KN so that I1 is smallest.
!
if ( i1 < i2 .and. i1 < i3 ) then
j = 3
else if ( i2 < i3 ) then
j = 2
isv = i1
i1 = i2
i2 = i3
i3 = isv
else
j = 1
isv = i1
i1 = i3
i3 = i2
i2 = isv
end if
!
! Test for KT < KN (triangle index not yet assigned).
!
if ( n1 < i1 ) then
cycle
end if
!
! Find KN, if it exists, by searching the triangle list in
! reverse order.
!
do kn = kt-1, 1, -1
if ( ltri(1,kn) == i1 .and. &
ltri(2,kn) == i2 .and. &
ltri(3,kn) == i3 ) then
go to 5
end if
end do
cycle
!
! Store KT as a neighbor of KN.
!
5 continue
ltri(j+3,kn) = kt
!
! Store KN as a neighbor of KT, and add a new arc KA.
!
6 continue
ltri(i+3,kt) = kn
if ( arcs ) then
ka = ka + 1
ltri(i+6,kt) = ka
if ( kn /= 0 ) then
ltri(j+6,kn) = ka
end if
end if
end do
!
! Bottom of loop on triangles.
!
8 continue
if ( lp2 /= lpln1 ) then
go to 1
end if
!9 continue
end do
nt = kt
ier = 0
return
end
subroutine trlist2 ( n, list, lptr, lend, nt, ltri, ier )
!*****************************************************************************80
!
!! TRLIST2 converts a triangulation data structure to a triangle list.
!
! Discussion:
!
! This subroutine converts a triangulation data structure
! from the linked list created by TRMESH to a triangle list.
!
! It is a version of TRLIST for the special case where the triangle
! list should only include the nodes that define each triangle.
!
! Modified:
!
! 21 July 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer N, the number of nodes in the triangulation.
! 3 <= N.
!
! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), linked
! list data structure defining the triangulation. Refer to TRMESH.
!
! Output, integer NT, the number of triangles in the
! triangulation unless IER /=0, in which case NT = 0. NT = 2N-NB-2 if
! NB >= 3 or 2N-4 if NB = 0, where NB is the number of boundary nodes.
!
! Output, integer LTRI(3,*). The second dimension of LTRI
! must be at least NT, where NT will be at most 2*N-4. The J-th column
! contains the vertex nodal indexes associated with triangle J for
! J = 1,...,NT. The vertices are ordered counterclockwise with the first
! vertex taken to be the one with smallest index. Thus, LTRI(2,J) and
! LTRI(3,J) are larger than LTRI(1,J) and index adjacent neighbors of node
! LTRI(1,J). The triangles are ordered on first (smallest) vertex indexes.
!
! Output, integer IER, error indicator.
! 0, if no errors were encountered.
! 1, if N is outside its valid range on input.
! 2, if the triangulation data structure (LIST,LPTR,LEND) is invalid.
! Note, however, that these arrays are not completely tested for validity.
!
! Local parameters:
!
! I,J = LTRI row indexes (1 to 3) associated with
! triangles KT and KN, respectively
! I1,I2,I3 = Nodal indexes of triangle KN
! ISV = Variable used to permute indexes I1,I2,I3
! KA = Arc index and number of currently stored arcs
! KN = Index of the triangle that shares arc I1-I2 with KT
! KT = Triangle index and number of currently stored triangles
! LP = LIST pointer
! LP2 = Pointer to N2 as a neighbor of N1
! LPL = Pointer to the last neighbor of I1
! LPLN1 = Pointer to the last neighbor of N1
! N1,N2,N3 = Nodal indexes of triangle KT
! NM2 = N-2
!
implicit none
integer n
integer i
integer i1
integer i2
integer i3
integer ier
integer isv
integer j
integer ka
integer kn
integer kt
integer lend(n)
integer list(6*(n-2))
integer lp
integer lp2
integer lpl
integer lpln1
integer lptr(6*(n-2))
integer ltri(3,*)
integer n1
integer n2
integer n3
integer nm2
integer nt
!
! Test for invalid input parameters.
!
if ( n < 3 ) then
nt = 0
ier = 1
return
end if
!
! Initialize parameters for loop on triangles KT = (N1,N2,
! N3), where N1 < N2 and N1 < N3.
!
! KA,KT = Numbers of currently stored arcs and triangles.
! NM2 = Upper bound on candidates for N1.
!
ka = 0
kt = 0
nm2 = n-2
!
! Loop on nodes N1.
!
do n1 = 1, nm2
!
! Loop on pairs of adjacent neighbors (N2,N3). LPLN1 points
! to the last neighbor of N1, and LP2 points to N2.
!
lpln1 = lend(n1)
lp2 = lpln1
1 continue
lp2 = lptr(lp2)
n2 = list(lp2)
lp = lptr(lp2)
n3 = abs ( list(lp) )
if ( n2 < n1 .or. n3 < n1 ) then
go to 8
end if
!
! Add a new triangle KT = (N1,N2,N3).
!
kt = kt + 1
ltri(1,kt) = n1
ltri(2,kt) = n2
ltri(3,kt) = n3
!
! Loop on triangle sides (I2,I1) with neighboring triangles
! KN = (I1,I2,I3).
!
do i = 1, 3
if ( i == 1 ) then
i1 = n3
i2 = n2
else if ( i == 2 ) then
i1 = n1
i2 = n3
else
i1 = n2
i2 = n1
end if
!
! Set I3 to the neighbor of I1 that follows I2 unless
! I2->I1 is a boundary arc.
!
lpl = lend(i1)
lp = lptr(lpl)
do
if ( list(lp) == i2 ) then
go to 3
end if
lp = lptr(lp)
if ( lp == lpl ) then
exit
end if
end do
!
! Invalid triangulation data structure: I1 is a neighbor of
! I2, but I2 is not a neighbor of I1.
!
if ( abs ( list(lp) ) /= i2 ) then
nt = 0
ier = 2
return
end if
!
! I2 is the last neighbor of I1. Bypass the search for a neighboring
! triangle if I2->I1 is a boundary arc.
!
kn = 0
if ( list(lp) < 0 ) then
go to 6
end if
!
! I2->I1 is not a boundary arc, and LP points to I2 as
! a neighbor of I1.
!
3 continue
lp = lptr(lp)
i3 = abs ( list(lp) )
!
! Find J such that LTRI(J,KN) = I3 (not used if KT < KN),
! and permute the vertex indexes of KN so that I1 is smallest.
!
if ( i1 < i2 .and. i1 < i3 ) then
j = 3
else if ( i2 < i3 ) then
j = 2
isv = i1
i1 = i2
i2 = i3
i3 = isv
else
j = 1
isv = i1
i1 = i3
i3 = i2
i2 = isv
end if
!
! Test for KT < KN (triangle index not yet assigned).
!
if ( n1 < i1 ) then
cycle
end if
!
! Find KN, if it exists, by searching the triangle list in
! reverse order.
!
do kn = kt-1, 1, -1
if ( ltri(1,kn) == i1 .and. &
ltri(2,kn) == i2 .and. &
ltri(3,kn) == i3 ) then
go to 5
end if
end do
cycle
5 continue
6 continue
end do
!
! Bottom of loop on triangles.
!
8 continue
if ( lp2 /= lpln1 ) then
go to 1
end if
!9 continue
end do
nt = kt
ier = 0
return
end
subroutine trlprt ( n, x, y, z, iflag, nrow, nt, ltri )
!*****************************************************************************80
!
!! TRLPRT prints a triangle list.
!
! Discussion:
!
! This subroutine prints the triangle list created by TRLIST
! and, optionally, the nodal coordinates
! (either latitude and longitude or Cartesian coordinates).
! The numbers of boundary nodes, triangles, and arcs are also printed.
!
! Modified:
!
! 06 June 2002
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer N, the number of nodes in the triangulation.
! 3 <= N <= 9999.
!
! Input, real ( kind = rk ) X(N), Y(N), Z(N), the coordinates of the nodes if
! IFLAG = 0, or (X and Y only) longitude and latitude, respectively,
! if 0 < IFLAG, or unused dummy parameters if IFLAG < 0.
!
! Input, integer IFLAG, nodal coordinate option indicator:
! = 0, if X, Y, and Z (assumed to contain Cartesian coordinates) are to
! be printed (to 6 decimal places).
! > 0, if only X and Y (assumed to contain longitude and latitude) are
! to be printed (to 6 decimal places).
! < 0, if only the adjacency lists are to be printed.
!
! Input, integer NROW, the number of rows (entries per triangle)
! reserved for the triangle list LTRI. The value must be 6 if only the
! vertex indexes and neighboring triangle indexes are stored, or 9
! if arc indexes are also stored.
!
! Input, integer NT, the number of triangles in the
! triangulation. 1 <= NT <= 9999.
!
! Input, integer LTRI(NROW,NT), the J-th column contains the
! vertex nodal indexes (first three rows), neighboring triangle indexes
! (second three rows), and, if NROW = 9, arc indexes (last three rows)
! associated with triangle J for J = 1,...,NT.
!
! Local parameters:
!
! I = DO-loop, nodal index, and row index for LTRI
! K = DO-loop and triangle index
! NA = Number of triangulation arcs
! NB = Number of boundary nodes
! NL = Number of lines printed on the current page
! NLMAX = Maximum number of print lines per page (except
! for the last page which may have two additional lines)
! NMAX = Maximum value of N and NT (4-digit format)
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer n
integer nrow
integer nt
integer i
integer iflag
integer k
integer ltri(nrow,nt)
integer na
integer nb
integer nl
integer, parameter :: nlmax = 58
integer, parameter :: nmax = 9999
real ( kind = rk ) x(n)
real ( kind = rk ) y(n)
real ( kind = rk ) z(n)
!
! Print a heading and test for invalid input.
!
write (*,100) n
nl = 3
if ( n < 3 .or. nmax < n .or. &
( nrow /= 6 .and. nrow /= 9) .or. &
nt < 1 .or. nmax < nt ) then
write (*,110) n, nrow, nt
return
end if
!
! Print X, Y, and Z.
!
if ( iflag == 0 ) then
write (*,101)
nl = 6
do i = 1, n
if ( nlmax <= nl ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
nl = 0
end if
write (*,103) i, x(i), y(i), z(i)
nl = nl + 1
end do
!
! Print X (longitude) and Y (latitude).
!
else if ( 0 < iflag ) then
write ( *, 102 )
nl = 6
do i = 1, n
if ( nlmax <= nl ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
nl = 0
end if
write (*,104) i, x(i), y(i)
nl = nl + 1
end do
end if
!
! Print the triangulation LTRI.
!
if ( nlmax / 2 < nl ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
nl = 0
end if
if ( nrow == 6 ) then
write (*,105)
else
write (*,106)
end if
nl = nl + 5
do k = 1, nt
if ( nlmax <= nl ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
nl = 0
end if
write (*,107) k, ltri(1:nrow,k)
nl = nl + 1
end do
!
! Print NB, NA, and NT (boundary nodes, arcs, and triangles).
!
nb = 2 * n - nt - 2
if ( nb < 3 ) then
nb = 0
na = 3 * n - 6
else
na = nt + n - 1
end if
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Number of boundary nodes NB = ', nb
write ( *, '(a,i8)' ) ' Number of arcs NA = ', na
write ( *, '(a,i8)' ) ' Number of triangles NT = ', nt
return
!
! Print formats:
!
100 format (///18x,'STRIPACK (TRLIST) output, n = ',i4)
101 format (//8x,'Node',10x,'X(node)',10x,'Y(node)',10x, &
'Z(node)'//)
102 format (//16x,'Node Longitude Latitude'//)
103 format (8x,i4,3e17.6)
104 format (16x,i4,2e17.6)
105 format (//' triangle',8x,'vertices',12x,'neighbors'/ &
4x,'kt',7x,'n1',5x,'n2',5x,'n3',4x,'kt1',4x, &
'kt2',4x,'kt3'/)
106 format (//'triangle',8x,'vertices',12x,'neighbors',14x,'arcs'/ &
4x,'kt n1 n2 n3 kt1',4x, &
'kt2 kt3 ka1 ka2 ka3'/)
107 format (2x,i4,2x,6(3x,i4),3(2x,i5))
110 format (//1x,10x,'Invalid parameter: N =',i5, &
', nrow =',i5,', nt =',i5,' ***')
end
subroutine trmesh ( n, x, y, z, list, lptr, lend, lnew, near, next, dist, ier )
!*****************************************************************************80
!
!! TRMESH creates a Delaunay triangulation on the unit sphere.
!
! Discussion:
!
! This subroutine creates a Delaunay triangulation of a
! set of N arbitrarily distributed points, referred to as
! nodes, on the surface of the unit sphere. The Delaunay
! triangulation is defined as a set of (spherical) triangles
! with the following five properties:
!
! 1) The triangle vertices are nodes.
! 2) No triangle contains a node other than its vertices.
! 3) The interiors of the triangles are pairwise disjoint.
! 4) The union of triangles is the convex hull of the set
! of nodes (the smallest convex set that contains
! the nodes). If the nodes are not contained in a
! single hemisphere, their convex hull is the
! entire sphere and there are no boundary nodes.
! Otherwise, there are at least three boundary nodes.
! 5) The interior of the circumcircle of each triangle
! contains no node.
!
! The first four properties define a triangulation, and the
! last property results in a triangulation which is as close
! as possible to equiangular in a certain sense and which is
! uniquely defined unless four or more nodes lie in a common
! plane. This property makes the triangulation well-suited
! for solving closest-point problems and for triangle-based
! interpolation.
!
! Provided the nodes are randomly ordered, the algorithm
! has expected time complexity O(N*log(N)) for most nodal
! distributions. Note, however, that the complexity may be
! as high as O(N**2) if, for example, the nodes are ordered
! on increasing latitude.
!
! Spherical coordinates (latitude and longitude) may be
! converted to Cartesian coordinates by Subroutine TRANS.
!
! The following is a list of the software package modules
! which a user may wish to call directly:
!
! ADDNOD - Updates the triangulation by appending a new node.
!
! AREAS - Returns the area of a spherical triangle.
!
! BNODES - Returns an array containing the indexes of the
! boundary nodes (if any) in counterclockwise
! order. Counts of boundary nodes, triangles,
! and arcs are also returned.
!
! CIRCUM - Returns the circumcenter of a spherical triangle.
!
! CRLIST - Returns the set of triangle circumcenters
! (Voronoi vertices) and circumradii associated
! with a triangulation.
!
! DELARC - Deletes a boundary arc from a triangulation.
!
! DELNOD - Updates the triangulation with a nodal deletion.
!
! EDGE - Forces an arbitrary pair of nodes to be connected
! by an arc in the triangulation.
!
! GETNP - Determines the ordered sequence of L closest nodes
! to a given node, along with the associated distances.
!
! INSIDE - Locates a point relative to a polygon on the
! surface of the sphere.
!
! INTRSC - Returns the point of intersection between a
! pair of great circle arcs.
!
! JRAND - Generates a uniformly distributed pseudo-random integer.
!
! LEFT - Locates a point relative to a great circle.
!
! NEARND - Returns the index of the nearest node to an
! arbitrary point, along with its squared
! distance.
!
! SCOORD - Converts a point from Cartesian coordinates to
! spherical coordinates.
!
! STORE - Forces a value to be stored in main memory so
! that the precision of floating point numbers
! in memory locations rather than registers is
! computed.
!
! TRANS - Transforms spherical coordinates into Cartesian
! coordinates on the unit sphere for input to
! Subroutine TRMESH.
!
! TRLIST - Converts the triangulation data structure to a
! triangle list more suitable for use in a finite
! element code.
!
! TRLPRT - Prints the triangle list created by TRLIST.
!
! TRMESH - Creates a Delaunay triangulation of a set of
! nodes.
!
! TRPLOT - Creates a level-2 Encapsulated Postscript (EPS)
! file containing a triangulation plot.
!
! TRPRNT - Prints the triangulation data structure and,
! optionally, the nodal coordinates.
!
! VRPLOT - Creates a level-2 Encapsulated Postscript (EPS)
! file containing a Voronoi diagram plot.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer N, the number of nodes in the triangulation.
! 3 <= N.
!
! Input, real ( kind = rk ) X(N), Y(N), Z(N), the coordinates of distinct
! nodes. (X(K),Y(K), Z(K)) is referred to as node K, and K is referred
! to as a nodal index. It is required that X(K)**2 + Y(K)**2 + Z(K)**2 = 1
! for all K. The first three nodes must not be collinear (lie on a
! common great circle).
!
! Output, integer LIST(6*(N-2)), nodal indexes which, along
! with LPTR, LEND, and LNEW, define the triangulation as a set of N
! adjacency lists; counterclockwise-ordered sequences of neighboring nodes
! such that the first and last neighbors of a boundary node are boundary
! nodes (the first neighbor of an interior node is arbitrary). In order to
! distinguish between interior and boundary nodes, the last neighbor of
! each boundary node is represented by the negative of its index.
!
! Output, integer LPTR(6*(N-2)), = Set of pointers (LIST
! indexes) in one-to-one correspondence with the elements of LIST.
! LIST(LPTR(I)) indexes the node which follows LIST(I) in cyclical
! counterclockwise order (the first neighbor follows the last neighbor).
!
! Output, integer LEND(N), pointers to adjacency lists.
! LEND(K) points to the last neighbor of node K. LIST(LEND(K)) < 0 if and
! only if K is a boundary node.
!
! Output, integer LNEW, pointer to the first empty location
! in LIST and LPTR (list length plus one). LIST, LPTR, LEND, and LNEW are
! not altered if IER < 0, and are incomplete if 0 < IER.
!
! Workspace, integer NEAR(N),
! used to efficiently determine the nearest triangulation node to each
! unprocessed node for use by ADDNOD.
!
! Workspace, integer NEXT(N),
! used to efficiently determine the nearest triangulation node to each
! unprocessed node for use by ADDNOD.
!
! Workspace, real ( kind = rk ) DIST(N),
! used to efficiently determine the nearest triangulation node to each
! unprocessed node for use by ADDNOD.
!
! Output, integer IER, error indicator:
! 0, if no errors were encountered.
! -1, if N < 3 on input.
! -2, if the first three nodes are collinear.
! L, if nodes L and M coincide for some L < M. The data structure
! represents a triangulation of nodes 1 to M-1 in this case.
!
! Local parameters:
!
! D = (Negative cosine of) distance from node K to node I
! D1,D2,D3 = Distances from node K to nodes 1, 2, and 3, respectively
! I,J = Nodal indexes
! I0 = Index of the node preceding I in a sequence of
! unprocessed nodes: I = NEXT(I0)
! K = Index of node to be added and DO-loop index: 3 < K
! LP = LIST index (pointer) of a neighbor of K
! LPL = Pointer to the last neighbor of K
! NEXTI = NEXT(I)
! NN = Local copy of N
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer n
real ( kind = rk ) d
real ( kind = rk ) d1
real ( kind = rk ) d2
real ( kind = rk ) d3
real ( kind = rk ) dist(n)
integer i
integer i0
integer ier
integer j
integer k
logical left
integer lend(n)
integer list(6*(n-2))
integer lnew
integer lp
integer lpl
integer lptr(6*(n-2))
integer near(n)
integer next(n)
integer nexti
integer nn
real ( kind = rk ) x(n)
real ( kind = rk ) y(n)
real ( kind = rk ) z(n)
nn = n
if ( nn < 3 ) then
ier = -1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRMESH - Fatal error!'
write ( *, '(a)' ) ' N < 3.'
stop
end if
!
! Store the first triangle in the linked list.
!
if ( .not. left (x(1),y(1),z(1),x(2),y(2),z(2), &
x(3),y(3),z(3) ) ) then
!
! The first triangle is (3,2,1) = (2,1,3) = (1,3,2).
!
list(1) = 3
lptr(1) = 2
list(2) = -2
lptr(2) = 1
lend(1) = 2
list(3) = 1
lptr(3) = 4
list(4) = -3
lptr(4) = 3
lend(2) = 4
list(5) = 2
lptr(5) = 6
list(6) = -1
lptr(6) = 5
lend(3) = 6
else if ( .not. left ( x(2),y(2),z(2),x(1),y(1),z(1),x(3),y(3),z(3) ) ) then
!
! The first triangle is (1,2,3): 3 Strictly Left 1->2,
! i.e., node 3 lies in the left hemisphere defined by arc 1->2.
!
list(1) = 2
lptr(1) = 2
list(2) = -3
lptr(2) = 1
lend(1) = 2
list(3) = 3
lptr(3) = 4
list(4) = -1
lptr(4) = 3
lend(2) = 4
list(5) = 1
lptr(5) = 6
list(6) = -2
lptr(6) = 5
lend(3) = 6
!
! The first three nodes are collinear.
!
else
ier = -2
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRMESH - Fatal error!'
write ( *, '(a)' ) ' The first 3 nodes are collinear.'
write ( *, '(a)' ) ' Try reordering the data.'
stop
end if
!
! Initialize LNEW and test for N = 3.
!
lnew = 7
if ( nn == 3 ) then
ier = 0
return
end if
!
! A nearest-node data structure (NEAR, NEXT, and DIST) is
! used to obtain an expected-time (N*log(N)) incremental
! algorithm by enabling constant search time for locating
! each new node in the triangulation.
!
! For each unprocessed node K, NEAR(K) is the index of the
! triangulation node closest to K (used as the starting
! point for the search in Subroutine TRFIND) and DIST(K)
! is an increasing function of the arc length (angular
! distance) between nodes K and NEAR(K): -Cos(a) for arc
! length a.
!
! Since it is necessary to efficiently find the subset of
! unprocessed nodes associated with each triangulation
! node J (those that have J as their NEAR entries), the
! subsets are stored in NEAR and NEXT as follows: for
! each node J in the triangulation, I = NEAR(J) is the
! first unprocessed node in J's set (with I = 0 if the
! set is empty), L = NEXT(I) (if 0 < I) is the second,
! NEXT(L) (if 0 < L) is the third, etc. The nodes in each
! set are initially ordered by increasing indexes (which
! maximizes efficiency) but that ordering is not main-
! tained as the data structure is updated.
!
! Initialize the data structure for the single triangle.
!
near(1) = 0
near(2) = 0
near(3) = 0
do k = nn, 4, -1
d1 = -( x(k) * x(1) + y(k) * y(1) + z(k) * z(1) )
d2 = -( x(k) * x(2) + y(k) * y(2) + z(k) * z(2) )
d3 = -( x(k) * x(3) + y(k) * y(3) + z(k) * z(3) )
if ( d1 <= d2 .and. d1 <= d3 ) then
near(k) = 1
dist(k) = d1
next(k) = near(1)
near(1) = k
else if ( d2 <= d1 .and. d2 <= d3 ) then
near(k) = 2
dist(k) = d2
next(k) = near(2)
near(2) = k
else
near(k) = 3
dist(k) = d3
next(k) = near(3)
near(3) = k
end if
end do
!
! Add the remaining nodes.
!
do k = 4, nn
call addnod ( near(k), k, x, y, z, list, lptr, lend, lnew, ier )
if ( ier /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRMESH - Fatal error!'
write ( *, '(a,i8)' ) ' ADDNOD returned error code IER = ', ier
stop
end if
!
! Remove K from the set of unprocessed nodes associated with NEAR(K).
!
i = near(k)
if ( near(i) == k ) then
near(i) = next(k)
else
i = near(i)
do
i0 = i
i = next(i0)
if ( i == k ) then
exit
end if
end do
next(i0) = next(k)
end if
near(k) = 0
!
! Loop on neighbors J of node K.
!
lpl = lend(k)
lp = lpl
3 continue
lp = lptr(lp)
j = abs ( list(lp) )
!
! Loop on elements I in the sequence of unprocessed nodes
! associated with J: K is a candidate for replacing J
! as the nearest triangulation node to I. The next value
! of I in the sequence, NEXT(I), must be saved before I
! is moved because it is altered by adding I to K's set.
!
i = near(j)
do
if ( i == 0 ) then
exit
end if
nexti = next(i)
!
! Test for the distance from I to K less than the distance
! from I to J.
!
d = - ( x(i) * x(k) + y(i) * y(k) + z(i) * z(k) )
if ( d < dist(i) ) then
!
! Replace J by K as the nearest triangulation node to I:
! update NEAR(I) and DIST(I), and remove I from J's set
! of unprocessed nodes and add it to K's set.
!
near(i) = k
dist(i) = d
if ( i == near(j) ) then
near(j) = nexti
else
next(i0) = nexti
end if
next(i) = near(k)
near(k) = i
else
i0 = i
end if
i = nexti
end do
!
! Bottom of loop on neighbors J.
!
!5 continue
if ( lp /= lpl ) then
go to 3
end if
!6 continue
end do
return
end
subroutine trplot ( lun, pltsiz, elat, elon, a, n, x, y, z, list, lptr, &
lend, title, numbr, ier )
!*****************************************************************************80
!
!! TRPLOT makes a PostScript image of a triangulation on a unit sphere.
!
! Discussion:
!
! This subroutine creates a level-2 Encapsulated Postscript (EPS)
! file containing a graphical display of a triangulation of a set of
! nodes on the unit sphere. The visible nodes are projected onto the
! plane that contains the origin and has normal defined by a
! user-specified eye-position. Projections of adjacent (visible) nodes
! are connected by line segments.
!
! The values in the data statements may be altered
! in order to modify various plotting options.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer LUN, the logical unit number in the range 0
! to 99. The unit should be opened with an appropriate
! file name before the call to this routine.
!
! Input, real ( kind = rk ) PLTSIZ, the plot size in inches. A circular
! window in the projection plane is mapped to a circular viewport with
! diameter equal to 0.88 * PLTSIZ (leaving room for labels outside the
! viewport). The viewport is centered on the 8.5 by 11 inch page, and
! its boundary is drawn. 1.0 <= PLTSIZ <= 8.5.
!
! Input, real ( kind = rk ) ELAT, ELON, the latitude and longitude
! (in degrees) of the center of projection E (the center of the plot).
! The projection plane is the plane that contains the origin and has
! E as unit normal. In a rotated coordinate system for which E is
! the north pole, the projection plane contains the equator, and only
! northern hemisphere nodes are visible (from the point at infinity in
! the direction E). These are projected orthogonally onto the
! projection plane (by zeroing the z-component in the rotated coordinate
! system). ELAT and ELON must be in the range -90 to 90 and -180 to
! 180, respectively.
!
! Input, real ( kind = rk ) A, the angular distance in degrees from E
! to the boundary of a circular window against which the triangulation
! is clipped. The projected window is a disk of radius R = Sin(A)
! centered at the origin, and only visible nodes whose projections are
! within distance R of the origin are included in the plot. Thus, if
! A = 90, the plot includes the entire hemisphere centered at E.
! 0 < A <= 90.
!
! Input, integer N, the number of nodes in the triangulation.
! 3 <= N.
!
! Input, real ( kind = rk ) X(N), Y(N), Z(N). the coordinates of the
! nodes (unit vectors).
!
! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), the
! data structure defining the triangulation, created by TRMESH.
!
! Input, character ( len = * ) TITLE, a string to be centered above the
! plot. The string must be enclosed in parentheses; i.e., the first and
! last characters must be '(' and ')', respectively, but these are not
! displayed. TITLE may have at most 80 characters including the parentheses.
!
! Input, logical NUMBR, option indicator: If NUMBR = TRUE, the
! nodal indexes are plotted next to the nodes.
!
! Output, integer IER, error indicator:
! 0, if no errors were encountered.
! 1, if LUN, PLTSIZ, or N is outside its valid range.
! 2, if ELAT, ELON, or A is outside its valid range.
! 3, if an error was encountered in writing to unit LUN.
!
! Local parameters:
!
! ANNOT = Logical variable with value TRUE iff the plot
! is to be annotated with the values of ELAT,
! ELON, and A
! CF = Conversion factor for degrees to radians
! CT = Cos(ELAT)
! EX,EY,EZ = Cartesian coordinates of the eye-position E
! FSIZN = Font size in points for labeling nodes with
! their indexes if NUMBR = TRUE
! FSIZT = Font size in points for the title (and
! annotation if ANNOT = TRUE)
! IPX1,IPY1 = X and y coordinates (in points) of the lower
! left corner of the bounding box or viewport box
! IPX2,IPY2 = X and y coordinates (in points) of the upper
! right corner of the bounding box or viewport box
! IR = Half the width (height) of the bounding box or
! viewport box in points -- viewport radius
! LP = LIST index (pointer)
! LPL = Pointer to the last neighbor of N0
! N0 = Index of a node whose incident arcs are to be drawn
! N1 = Neighbor of N0
! R11...R23 = Components of the first two rows of a rotation
! that maps E to the north pole (0,0,1)
! SF = Scale factor for mapping world coordinates
! (window coordinates in [-WR,WR] X [-WR,WR])
! to viewport coordinates in [IPX1,IPX2] X [IPY1,IPY2]
! T = Temporary variable
! TX,TY = Translation vector for mapping world coordi-
! nates to viewport coordinates
! WR = Window radius r = Sin(A)
! WRS = WR**2
! X0,Y0,Z0 = Coordinates of N0 in the rotated coordinate
! system or label location (X0,Y0)
! X1,Y1,Z1 = Coordinates of N1 in the rotated coordinate
! system or intersection of edge N0-N1 with
! the equator (in the rotated coordinate system)
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer n
real ( kind = rk ) a
logical, parameter :: annot = .true.
real ( kind = rk ) cf
real ( kind = rk ) ct
real ( kind = rk ) elat
real ( kind = rk ) elon
real ( kind = rk ) ex
real ( kind = rk ) ey
real ( kind = rk ) ez
real ( kind = rk ), parameter :: fsizn = 10.0D+00
real ( kind = rk ), parameter :: fsizt = 16.0D+00
integer ier
integer ipx1
integer ipx2
integer ipy1
integer ipy2
integer ir
integer lend(n)
integer list(6*(n-2))
integer lp
integer lpl
integer lptr(6*(n-2))
integer lun
integer n0
integer n1
logical numbr
real ( kind = rk ) pltsiz
real ( kind = rk ) r11
real ( kind = rk ) r12
real ( kind = rk ) r21
real ( kind = rk ) r22
real ( kind = rk ) r23
real ( kind = rk ) sf
real ( kind = rk ) t
character ( len = * ) title
real ( kind = rk ) tx
real ( kind = rk ) ty
real ( kind = rk ) wr
real ( kind = rk ) wrs
real ( kind = rk ) x(n)
real ( kind = rk ) x0
real ( kind = rk ) x1
real ( kind = rk ) y(n)
real ( kind = rk ) y0
real ( kind = rk ) y1
real ( kind = rk ) z(n)
real ( kind = rk ) z0
real ( kind = rk ) z1
ier = 0
!
! Test for invalid parameters.
!
if ( lun < 0 ) then
ier = 1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPLOT - Fatal error!'
write ( *, '(a)' ) ' LUN < 0!'
return
end if
if ( 99 < lun ) then
ier = 1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPLOT - Fatal error!'
write ( *, '(a)' ) ' 99 < LUN!'
return
end if
if ( pltsiz < 1.0D+00 ) then
ier = 1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPLOT - Fatal error!'
write ( *, '(a)' ) ' PLTSIZ < 1.0!'
return
else if ( 8.5D+00 < pltsiz ) then
ier = 1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPLOT - Fatal error!'
write ( *, '(a)' ) ' 8.5 < PLTSIZ!'
return
else if ( n < 3 ) then
ier = 1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPLOT - Fatal error!'
write ( *, '(a)' ) ' N < 3!'
return
end if
if ( 90.0D+00 < abs ( elat ) ) then
ier = 2
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPLOT - Fatal error!'
write ( *, '(a)' ) ' 90 < |ELAT|!'
return
else if ( 180.0D+00 < abs ( elon ) ) then
ier = 2
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPLOT - Fatal error!'
write ( *, '(a)' ) ' 180 < |ELON|!'
return
else if ( 90.0D+00 < a ) then
ier = 2
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPLOT - Fatal error!'
write ( *, '(a)' ) ' 90.0 < A!'
return
end if
!
! Compute a conversion factor CF for degrees to radians.
!
cf = atan ( 1.0D+00 ) / 45.0D+00
!
! Compute the window radius WR.
!
wr = sin ( cf * a )
wrs = wr * wr
!
! Compute the lower left (IPX1,IPY1) and upper right
! (IPX2,IPY2) corner coordinates of the bounding box.
! The coordinates, specified in default user space units
! (points, at 72 points/inch with origin at the lower
! left corner of the page), are chosen to preserve the
! square aspect ratio, and to center the plot on the 8.5
! by 11 inch page. The center of the page is (306,396),
! and IR = PLTSIZ/2 in points.
!
ir = nint ( 36.0D+00 * pltsiz )
ipx1 = 306 - ir
ipx2 = 306 + ir
ipy1 = 396 - ir
ipy2 = 396 + ir
!
! Output header comments.
!
write ( lun, '(a)' ) '%!ps-adobe-3.0 epsf-3.0'
write ( lun, '(a,4i4)' ) '%%BoundingBox:', ipx1, ipy1, ipx2, ipy2
write ( lun, '(a)' ) '%%title: Triangulation'
write ( lun, '(a)' ) '%%creator: STRIPACK.F90'
write ( lun, '(a)' ) '%%endcomments'
!
! Set (IPX1,IPY1) and (IPX2,IPY2) to the corner coordinates
! of a viewport box obtained by shrinking the bounding box
! by 12% in each dimension.
!
ir = nint ( 0.88D+00 * real ( ir, kind = rk ) )
ipx1 = 306 - ir
ipx2 = 306 + ir
ipy1 = 396 - ir
ipy2 = 396 + ir
!
! Set the line thickness to 2 points, and draw the
! viewport boundary.
!
t = 2.0D+00
write ( lun, '(f12.6,a)' ) t, ' setlinewidth'
write ( lun, '(a,i3,a)' ) '306 396 ', ir, ' 0 360 arc'
write ( lun, '(a)' ) 'stroke'
!
! Set up an affine mapping from the window box [-WR,WR] X
! [-WR,WR] to the viewport box.
!
sf = real ( ir, kind = rk ) / wr
tx = ipx1 + sf * wr
ty = ipy1 + sf * wr
write ( lun, '(2f12.6,a)' ) tx, ty, ' translate'
write ( lun, '(2f12.6,a)' ) sf, sf, ' scale'
!
! The line thickness must be changed to reflect the new
! scaling which is applied to all subsequent output.
! Set it to 1.0 point.
!
t = 1.0D+00 / sf
write ( lun, '(f12.6,a)' ) t, ' setlinewidth'
!
! Save the current graphics state, and set the clip path to
! the boundary of the window.
!
write ( lun, '(a)' ) 'gsave'
write ( lun, '(a,f12.6,a)' ) '0 0 ', wr, ' 0 360 arc'
write ( lun, '(a)' ) 'clip newpath'
!
! Compute the Cartesian coordinates of E and the components
! of a rotation R which maps E to the north pole (0,0,1).
! R is taken to be a rotation about the z-axis (into the
! yz-plane) followed by a rotation about the x-axis chosen
! so that the view-up direction is (0,0,1), or (-1,0,0) if
! E is the north or south pole.
!
! ( R11 R12 0 )
! R = ( R21 R22 R23 )
! ( EX EY EZ )
!
t = cf * elon
ct = cos ( cf * elat )
ex = ct * cos ( t )
ey = ct * sin ( t )
ez = sin ( cf * elat )
if ( ct /= 0.0D+00 ) then
r11 = -ey / ct
r12 = ex / ct
else
r11 = 0.0D+00
r12 = 1.0D+00
end if
r21 = -ez * r12
r22 = ez * r11
r23 = ct
!
! Loop on visible nodes N0 that project to points (X0,Y0) in the window.
!
do n0 = 1, n
z0 = ex * x(n0) + ey * y(n0) + ez * z(n0)
if ( z0 < 0.0D+00 ) then
cycle
end if
x0 = r11 * x(n0) + r12 * y(n0)
y0 = r21 * x(n0) + r22 * y(n0) + r23 * z(n0)
if ( wrs < x0 * x0 + y0 * y0 ) then
cycle
end if
lpl = lend(n0)
lp = lpl
!
! Loop on neighbors N1 of N0. LPL points to the last
! neighbor of N0. Copy the components of N1 into P.
!
do
lp = lptr(lp)
n1 = abs ( list(lp) )
x1 = r11 * x(n1) + r12 * y(n1)
y1 = r21 * x(n1) + r22 * y(n1) + r23 * z(n1)
z1 = ex * x(n1) + ey * y(n1) + ez * z(n1)
!
! N1 is a 'southern hemisphere' point. Move it to the
! intersection of edge N0-N1 with the equator so that
! the edge is clipped properly. Z1 is implicitly set
! to 0.
!
if ( z1 < 0.0D+00 ) then
x1 = z0 * x1 - z1 * x0
y1 = z0 * y1 - z1 * y0
t = sqrt ( x1 * x1 + y1 * y1 )
x1 = x1 / t
y1 = y1 / t
end if
!
! If node N1 is in the window and N1 < N0, bypass edge
! N0->N1 (since edge N1->N0 has already been drawn).
!
! Add the edge to the path.
!
if ( z1 < 0.0D+00 .or. &
wrs < x1 * x1 + y1 * y1 .or. &
n0 <= n1 ) then
write ( lun, '(2f12.6,a,2f12.6,a)' ) &
x0, y0, ' moveto', x1, y1, ' lineto'
end if
if ( lp == lpl ) then
exit
end if
end do
end do
!
! Paint the path and restore the saved graphics state (with
! no clip path).
!
write ( lun, '(a)' ) 'stroke'
write ( lun, '(a)' ) 'grestore'
if ( numbr ) then
!
! Nodes in the window are to be labeled with their indexes.
! Convert FSIZN from points to world coordinates, and
! output the commands to select a font and scale it.
!
t = fsizn / sf
write ( lun, '(a)' ) '/Helvetica findfont'
write ( lun, '(f12.6,a)' ) t, ' scalefont setfont'
!
! Loop on visible nodes N0 that project to points (X0,Y0) in the window.
!
do n0 = 1, n
if ( ex * x(n0) + ey * y(n0) + ez * z(n0) < 0.0D+00 ) then
cycle
end if
x0 = r11 * x(n0) + r12 * y(n0)
y0 = r21 * x(n0) + r22 * y(n0) + r23 * z(n0)
if ( wrs < x0 * x0 + y0 * y0 ) then
cycle
end if
!
! Move to (X0,Y0) and draw the label N0. The first char-
! acter will will have its lower left corner about one
! character width to the right of the nodal position.
!
write ( lun, '(2f12.6,a)' ) x0, y0, ' moveto'
write ( lun, '(a,i3,a)' ) '(', n0, ') show'
end do
end if
!
! Convert FSIZT from points to world coordinates, and output
! the commands to select a font and scale it.
!
t = fsizt / sf
write ( lun, '(a)' ) '/Helvetica findfont'
write ( lun, '(f12.6,a)' ) t, ' scalefont setfont'
!
! Display TITLE centered above the plot:
!
y0 = wr + 3.0D+00 * t
write ( lun, '(a)' ) title
write ( lun, '(a,f12.6,a)' ) ' stringwidth pop 2 div neg ', y0, ' moveto'
write ( lun, '(a)' ) title
write ( lun, '(a)' ) ' show'
!
! Display the window center and radius below the plot.
!
if ( annot ) then
x0 = -wr
y0 = -wr - 50.0D+00 / sf
write ( lun, '(2f12.6,a)' ) x0, y0, ' moveto'
write ( lun, '(a,f7.2,a,f8.2,a)' ) '(Window center: Latitude = ', elat, &
', Longitude = ', elon , ') show'
y0 = y0 - 2.0D+00 * t
write ( lun, '(2f12.6,a)' ) x0, y0, ' moveto'
write ( lun, '(a,f5.2,a)' ) '(Angular extent = ', a, ') show'
end if
!
! Paint the path and output the showpage command and
! end-of-file indicator.
!
write ( lun, '(a)' ) 'stroke'
write ( lun, '(a)' ) 'showpage'
write ( lun, '(a)' ) '%%eof'
ier = 0
return
end
subroutine trprnt ( n, x, y, z, iflag, list, lptr, lend )
!*****************************************************************************80
!
!! TRPRNT prints the triangulation adjacency lists.
!
! Discussion:
!
! This subroutine prints the triangulation adjacency lists
! created by TRMESH and, optionally, the nodal
! coordinates (either latitude and longitude or Cartesian
! coordinates) on logical unit LOUT. The list of neighbors
! of a boundary node is followed by index 0. The numbers of
! boundary nodes, triangles, and arcs are also printed.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer N = Number of nodes in the triangulation.
! 3 <= N and N <= 9999.
!
! Input, real ( kind = rk ) X(N), Y(N), Z(N), the Cartesian coordinates of
! the nodes if IFLAG = 0, or (X and Y only) containing longitude and
! latitude, respectively, if 0 < IFLAG, or unused dummy parameters if
! IFLAG < 0.
!
! Input, integer IFLAG = Nodal coordinate option indicator:
! = 0 if X, Y, and Z (assumed to contain Cartesian coordinates) are to be
! printed (to 6 decimal places).
! > 0 if only X and Y (assumed to contain longitude and latitude) are
! to be printed (to 6 decimal places).
! < 0 if only the adjacency lists are to be printed.
!
! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), the
! data structure defining the triangulation. Refer to TRMESH.
!
! Local parameters:
!
! I = NABOR index (1 to K)
! INC = Increment for NL associated with an adjacency list
! K = Counter and number of neighbors of NODE
! LP = LIST pointer of a neighbor of NODE
! LPL = Pointer to the last neighbor of NODE
! NA = Number of arcs in the triangulation
! NABOR = Array containing the adjacency list associated
! with NODE, with zero appended if NODE is a boundary node
! NB = Number of boundary nodes encountered
! ND = Index of a neighbor of NODE (or negative index)
! NL = Number of lines that have been printed on the current page
! NLMAX = Maximum number of print lines per page (except
! for the last page which may have two additional lines)
! NMAX = Upper bound on N (allows 4-digit indexes)
! NODE = Index of a node and DO-loop index (1 to N)
! NN = Local copy of N
! NT = Number of triangles in the triangulation
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer n
integer iflag
integer inc
integer k
integer lend(n)
integer list(6*(n-2))
integer lp
integer lpl
integer lptr(6*(n-2))
integer na
integer nabor(400)
integer nb
integer nd
integer nl
integer, parameter :: nlmax = 58
integer, parameter :: nmax = 9999
integer nn
integer node
integer nt
real ( kind = rk ) x(n)
real ( kind = rk ) y(n)
real ( kind = rk ) z(n)
nn = n
!
! Print a heading and test the range of N.
!
write (*,100) nn
if ( nn < 3 .or. nmax < nn ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TRPRNT - Fatal error!'
write ( *, '(a)' ) ' N is outside its valid range.'
return
end if
!
! Initialize NL (the number of lines printed on the current
! page) and NB (the number of boundary nodes encountered).
!
nl = 6
nb = 0
!
! Print LIST only. K is the number of neighbors of NODE
! that have been stored in NABOR.
!
if ( iflag < 0 ) then
write (*,101)
do node = 1, nn
lpl = lend(node)
lp = lpl
k = 0
do
k = k + 1
lp = lptr(lp)
nd = list(lp)
nabor(k) = nd
if ( lp == lpl ) then
exit
end if
end do
!
! NODE is a boundary node. Correct the sign of the last
! neighbor, add 0 to the end of the list, and increment NB.
!
if ( nd <= 0 ) then
nabor(k) = -nd
k = k + 1
nabor(k) = 0
nb = nb + 1
end if
!
! Increment NL and print the list of neighbors.
!
inc = ( k - 1 ) / 14 + 2
nl = nl + inc
if ( nlmax < nl ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
nl = inc
end if
write (*,104) node, nabor(1:k)
if ( k /= 14 ) then
write ( *, '(a)' ) ' '
end if
end do
else if ( 0 < iflag ) then
!
! Print X (longitude), Y (latitude), and LIST.
!
write (*,102)
do node = 1, nn
lpl = lend(node)
lp = lpl
k = 0
do
k = k + 1
lp = lptr(lp)
nd = list(lp)
nabor(k) = nd
if ( lp == lpl ) then
exit
end if
end do
if ( nd <= 0 ) then
!
! NODE is a boundary node.
!
nabor(k) = -nd
k = k + 1
nabor(k) = 0
nb = nb + 1
end if
!
! Increment NL and print X, Y, and NABOR.
!
inc = ( k - 1 ) / 8 + 2
nl = nl + inc
if ( nlmax < nl ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
nl = inc
end if
write (*,105) node, x(node), y(node), nabor(1:k)
if ( k /= 8 ) then
write ( *, '(a)' ) ' '
end if
end do
else
!
! Print X, Y, Z, and LIST.
!
write (*,103)
do node = 1, nn
lpl = lend(node)
lp = lpl
k = 0
do
k = k + 1
lp = lptr(lp)
nd = list(lp)
nabor(k) = nd
if ( lp == lpl ) then
exit
end if
end do
!
! NODE is a boundary node.
!
if ( nd <= 0 ) then
nabor(k) = -nd
k = k + 1
nabor(k) = 0
nb = nb + 1
end if
!
! Increment NL and print X, Y, Z, and NABOR.
!
inc = ( k - 1 ) / 5 + 2
nl = nl + inc
if ( nlmax < nl ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' '
nl = inc
end if
write (*,106) node, x(node), y(node), z(node), nabor(1:k)
if ( k /= 5 ) then
write ( *, '(a)' ) ' '
end if
end do
end if
!
! Print NB, NA, and NT (boundary nodes, arcs, and triangles).
!
if ( nb /= 0 ) then
na = 3 * nn - nb - 3
nt = 2 * nn - nb - 2
else
na = 3 * nn - 6
nt = 2 * nn - 4
end if
write ( *, '(a)' ) ' '
write ( *, '(a,i8,a)' ) ' NB = ', nb, ' boundary arcs.'
write ( *, '(a,i8,a)' ) ' NA = ', na, ' arcs.'
write ( *, '(a,i8,a)' ) ' NT = ', nt, ' triangles.'
return
!
! Print formats:
!
100 format (///15x,'STRIPACK triangulation data ', &
'structure, n = ',i5//)
101 format (' node',31x,'neighbors of node'//)
102 format (' Node Longitude Latitude', &
18x,'neighbors of node'//)
103 format (' node x(node) y(node)',8x, &
'z(node)',11x,'neighbors of node'//)
104 format (i5,4x,14i5/(1x,8x,14i5))
105 format (i5,2e15.6,4x,8i5/(1x,38x,8i5))
106 format (i5,3e15.6,4x,5i5/(1x,53x,5i5))
end
subroutine voronoi_poly_count ( n, lend, lptr, listc )
!*****************************************************************************80
!
!! VORONOI_POLY_COUNT counts the polygons of each size in the Voronoi diagram.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 06 June 2002
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the number of Voronoi polygons.
!
! Input, integer LEND(N), some kind of pointer.
!
! Input, integer LPTR(6*(N-2)), some other kind of pointer.
!
! Input, integer LISTC(6*(N-2)), some other kind of pointer.
!
implicit none
integer n
integer, parameter :: side_max = 20
integer count(side_max)
integer edges
integer i
integer kv
integer lend(n)
integer listc(6*(n-2))
integer lp
integer lpl
integer lptr(6*(n-2))
integer n0
integer sides
integer vertices
count(1:side_max) = 0
edges = 0
vertices = 0
do n0 = 1, n
lpl = lend(n0)
lp = lpl
sides = 0
do
lp = lptr(lp)
kv = listc(lp)
vertices = max ( vertices, kv )
sides = sides + 1
edges = edges + 1
if ( lp == lpl ) then
exit
end if
end do
if ( 0 < sides .and. sides < side_max ) then
count(sides) = count(sides) + 1
else
count(side_max) = count(side_max) + 1
end if
end do
edges = edges / 2
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VORONOI_POLY_COUNT:'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Check Euler''s formula:'
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Faces = ', n
write ( *, '(a,i8)' ) ' Vertices = ', vertices
write ( *, '(a,i8)' ) ' Edges = ', edges
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' F+V-E-2 = ', n + vertices - edges - 2
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Count the polygons having a given number of sides.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Sides Number'
write ( *, '(a)' ) ' '
do i = 1, side_max - 1
if ( count(i) /= 0 ) then
write ( *, '(2x,i8,2x,i8)' ) i, count(i)
end if
end do
if ( count(side_max) /= 0 ) then
write ( *, '(2x,i8,2x,i8)' ) side_max, count(side_max)
end if
return
end
subroutine vrplot ( lun, pltsiz, elat, elon, a, n, x, y, z, nt, listc, lptr, &
lend, xc, yc, zc, title, numbr, ier )
!*****************************************************************************80
!
!! VRPLOT makes a PostScript image of a Voronoi diagram on the unit sphere.
!
! Discussion:
!
! This subroutine creates a level-2 Encapsulated Postscript
! (EPS) file containing a graphical depiction of a
! Voronoi diagram of a set of nodes on the unit sphere.
! The visible vertices are projected onto the plane that
! contains the origin and has normal defined by a user-
! specified eye-position. Projections of adjacent (visible)
! Voronoi vertices are connected by line segments.
!
! The parameters defining the Voronoi diagram may be computed by
! subroutine CRLIST.
!
! Modified:
!
! 16 June 2007
!
! Author:
!
! Robert Renka
!
! Reference:
!
! Robert Renka,
! Algorithm 772: STRIPACK,
! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
! ACM Transactions on Mathematical Software,
! Volume 23, Number 3, September 1997, pages 416-434.
!
! Parameters:
!
! Input, integer LUN, the logical unit number in the range 0
! to 99. The unit should be opened with an appropriate
! file name before the call to this routine.
!
! Input, real ( kind = rk ) PLTSIZ, the plot size in inches. A circular
! window in the projection plane is mapped to a circular viewport with
! diameter equal to .88*PLTSIZ (leaving room for labels outside the
! viewport). The viewport is centered on the 8.5 by 11 inch page, and its
! boundary is drawn. 1.0 <= PLTSIZ <= 8.5.
!
! Input, real ( kind = rk ) ELAT, ELON, the latitude and longitude (in
! degrees) of the center of projection E (the center of the plot). The
! projection plane is the plane that contains the origin and has E as unit
! normal. In a rotated coordinate system for which E is the north pole, the
! projection plane contains the equator, and only northern hemisphere
! points are visible (from the point at infinity in the direction E).
! These are projected orthogonally onto the projection plane (by zeroing
! the z-component in the rotated coordinate system). ELAT and ELON must
! be in the range -90 to 90 and -180 to 180, respectively.
!
! Input, real ( kind = rk ) A, the angular distance in degrees from E to the
! boundary of a circular window against which the Voronoi diagram is clipped.
! The projected window is a disk of radius R = Sin(A) centered at the
! origin, and only visible vertices whose projections are within distance
! R of the origin are included in the plot. Thus, if A = 90, the plot
! includes the entire hemisphere centered at E. 0 < A <= 90.
!
! Input, integer N, the number of nodes (Voronoi centers) and
! Voronoi regions. 3 <= N.
!
! Input, real ( kind = rk ) X(N), Y(N), Z(N), the coordinates of the nodes
! (unit vectors).
!
! Input, integer NT, the number of Voronoi region vertices
! (triangles, including those in the extended triangulation if the number
! of boundary nodes NB is nonzero): NT = 2*N-4.
!
! Input, integer LISTC(3*NT), containing triangle indexes
! (indexes to XC, YC, and ZC) stored in 1-1 correspondence with LIST/LPTR
! entries (or entries that would be stored in LIST for the extended
! triangulation): the index of triangle (N1,N2,N3) is stored in LISTC(K),
! LISTC(L), and LISTC(M), where LIST(K), LIST(L), and LIST(M) are the
! indexes of N2 as a neighbor of N1, N3 as a neighbor of N2, and N1 as a
! neighbor of N3. The Voronoi region associated with a node is defined by
! the CCW-ordered sequence of circumcenters in one-to-one correspondence
! with its adjacency list (in the extended triangulation).
!
! Input, integer LPTR(3*NT), where NT = 2*N-4, containing a
! set of pointers (LISTC indexes) in one-to-one correspondence with the
! elements of LISTC. LISTC(LPTR(I)) indexes the triangle which follows
! LISTC(I) in cyclical counterclockwise order (the first neighbor follows
! the last neighbor).
!
! Input, integer LEND(N), a set of pointers to triangle lists.
! LP = LEND(K) points to a triangle (indexed by LISTC(LP)) containing node
! K for K = 1 to N.
!
! Input, real ( kind = rk ) XC(NT), YC(NT), ZC(NT), the coordinates of the
! triangle circumcenters (Voronoi vertices).
! XC(I)**2 + YC(I)**2 + ZC(I)**2 = 1.
!
! Input, character ( len = * ) TITLE, a string to be centered above the plot.
! The string must be enclosed in parentheses; i.e., the first and last
! characters must be '(' and ')', respectively, but these are not
! displayed. TITLE may have at most 80 characters including the parentheses.
!
! Input, logical NUMBR, option indicator: If NUMBR = TRUE, the nodal
! indexes are plotted at the Voronoi region centers.
!
! Output, integer IER = Error indicator:
! 0, if no errors were encountered.
! 1, if LUN, PLTSIZ, N, or NT is outside its valid range.
! 2, if ELAT, ELON, or A is outside its valid range.
! 3, if an error was encountered in writing to unit LUN.
!
! Local parameters:
!
! ANNOT = Logical variable with value TRUE iff the plot
! is to be annotated with the values of ELAT, ELON, and A
! CF = Conversion factor for degrees to radians
! CT = Cos(ELAT)
! EX,EY,EZ = Cartesian coordinates of the eye-position E
! FSIZN = Font size in points for labeling nodes with
! their indexes if NUMBR = TRUE
! FSIZT = Font size in points for the title (and
! annotation if ANNOT = TRUE)
! IN1,IN2 = Logical variables with value TRUE iff the
! projections of vertices KV1 and KV2, respec-
! tively, are inside the window
! IPX1,IPY1 = X and y coordinates (in points) of the lower
! left corner of the bounding box or viewport box
! IPX2,IPY2 = X and y coordinates (in points) of the upper
! right corner of the bounding box or viewport box
! IR = Half the width (height) of the bounding box or
! viewport box in points -- viewport radius
! KV1,KV2 = Endpoint indexes of a Voronoi edge
! LP = LIST index (pointer)
! LPL = Pointer to the last neighbor of N0
! N0 = Index of a node
! R11...R23 = Components of the first two rows of a rotation
! that maps E to the north pole (0,0,1)
! SF = Scale factor for mapping world coordinates
! (window coordinates in [-WR,WR] X [-WR,WR])
! to viewport coordinates in [IPX1,IPX2] X [IPY1,IPY2]
! T = Temporary variable
! TX,TY = Translation vector for mapping world coordi-
! nates to viewport coordinates
! WR = Window radius r = Sin(A)
! WRS = WR**2
! X0,Y0 = Projection plane coordinates of node N0 or label location
! X1,Y1,Z1 = Coordinates of vertex KV1 in the rotated coordinate system
! X2,Y2,Z2 = Coordinates of vertex KV2 in the rotated
! coordinate system or intersection of edge
! KV1-KV2 with the equator (in the rotated coordinate system)
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
integer n
integer nt
real ( kind = rk ) a
logical, parameter :: annot = .true.
real ( kind = rk ) cf
real ( kind = rk ) ct
real ( kind = rk ) elat
real ( kind = rk ) elon
real ( kind = rk ) ex
real ( kind = rk ) ey
real ( kind = rk ) ez
real ( kind = rk ), parameter :: fsizn = 10.0D+00
real ( kind = rk ), parameter :: fsizt = 16.0D+00
integer ier
logical in1
logical in2
integer ipx1
integer ipx2
integer ipy1
integer ipy2
integer ir
integer kv1
integer kv2
integer lend(n)
integer listc(3*nt)
integer lp
integer lpl
integer lptr(6*(n-2))
integer lun
integer n0
logical numbr
real ( kind = rk ) pltsiz
real ( kind = rk ) r11
real ( kind = rk ) r12
real ( kind = rk ) r21
real ( kind = rk ) r22
real ( kind = rk ) r23
real ( kind = rk ) sf
real ( kind = rk ) t
character ( len = * ) title
real ( kind = rk ) tx
real ( kind = rk ) ty
real ( kind = rk ) wr
real ( kind = rk ) wrs
real ( kind = rk ) x(n)
real ( kind = rk ) x0
real ( kind = rk ) x1
real ( kind = rk ) x2
real ( kind = rk ) xc(nt)
real ( kind = rk ) y(n)
real ( kind = rk ) y0
real ( kind = rk ) y1
real ( kind = rk ) y2
real ( kind = rk ) yc(nt)
real ( kind = rk ) z(n)
real ( kind = rk ) z1
real ( kind = rk ) z2
real ( kind = rk ) zc(nt)
ier = 0
!
! Test for invalid parameters.
!
if ( lun < 0 ) then
ier = 1
return
end if
if ( 99 < lun ) then
ier = 1
return
end if
if ( pltsiz < 1.0D+00 .or. 8.5D+00 < pltsiz .or. &
n < 3 .or. nt /= 2*n-4) then
ier = 1
return
end if
if ( 90.0D+00 < abs ( elat ) .or. &
180.0D+00 < abs ( elon ) .or. &
90.0D+00 < a ) then
ier = 2
return
end if
!
! Compute a conversion factor CF for degrees to radians.
!
cf = atan ( 1.0D+00 ) / 45.0D+00
!
! Compute the window radius WR.
!
wr = sin ( cf * a )
wrs = wr * wr
!
! Compute the lower left (IPX1,IPY1) and upper right
! (IPX2,IPY2) corner coordinates of the bounding box.
! The coordinates, specified in default user space units
! (points, at 72 points/inch with origin at the lower
! left corner of the page), are chosen to preserve the
! square aspect ratio, and to center the plot on the 8.5
! by 11 inch page. The center of the page is (306,396),
! and IR = PLTSIZ/2 in points.
!
ir = nint ( 36.0D+00 * pltsiz )
ipx1 = 306 - ir
ipx2 = 306 + ir
ipy1 = 396 - ir
ipy2 = 396 + ir
!
! Output header comments.
!
write ( lun, '(a)' ) '%!ps-adobe-3.0 epsf-3.0'
write ( lun, '(a,4i4)' ) '%%BoundingBox: ', ipx1, ipy1, ipx2, ipy2
write ( lun, '(a)' ) '%%title: Voronoi diagram'
write ( lun, '(a)' ) '%%creator: STRIPACK.F90'
write ( lun, '(a)' ) '%%endcomments'
!
! Set (IPX1,IPY1) and (IPX2,IPY2) to the corner coordinates
! of a viewport box obtained by shrinking the bounding box
! by 12% in each dimension.
!
ir = nint ( 0.88D+00 * real ( ir, kind = rk ) )
ipx1 = 306 - ir
ipx2 = 306 + ir
ipy1 = 396 - ir
ipy2 = 396 + ir
!
! Set the line thickness to 2 points, and draw the viewport boundary.
!
t = 2.0D+00
write ( lun, '(f12.6,a)' ) t, ' setlinewidth'
write ( lun, '(a,i3,a)' ) '306 396 ', ir, ' 0 360 arc'
write ( lun, '(a)' ) 'stroke'
!
! Set up an affine mapping from the window box [-WR,WR] X
! [-WR,WR] to the viewport box.
!
sf = real ( ir, kind = rk ) / wr
tx = ipx1 + sf * wr
ty = ipy1 + sf * wr
write ( lun, '(2f12.6,a)' ) tx, ty, ' translate'
write ( lun, '(2f12.6,a)' ) sf, sf, ' scale'
!
! The line thickness must be changed to reflect the new
! scaling which is applied to all subsequent output.
! Set it to 1.0 point.
!
t = 1.0D+00 / sf
write ( lun, '(f12.6,a)' ) t, ' setlinewidth'
!
! Save the current graphics state, and set the clip path to
! the boundary of the window.
!
write ( lun, '(a)' ) 'gsave'
write ( lun, '(a,f12.6,a)' ) '0 0 ', wr, ' 0 360 arc'
write ( lun, '(a)' ) 'clip newpath'
!
! Compute the Cartesian coordinates of E and the components
! of a rotation R which maps E to the north pole (0,0,1).
! R is taken to be a rotation about the z-axis (into the
! yz-plane) followed by a rotation about the x-axis chosen
! so that the view-up direction is (0,0,1), or (-1,0,0) if
! E is the north or south pole.
!
! ( R11 R12 0 )
! R = ( R21 R22 R23 )
! ( EX EY EZ )
!
t = cf * elon
ct = cos ( cf * elat )
ex = ct * cos ( t )
ey = ct * sin ( t )
ez = sin ( cf * elat )
if ( ct /= 0.0D+00 ) then
r11 = -ey / ct
r12 = ex / ct
else
r11 = 0.0D+00
r12 = 1.0D+00
end if
r21 = -ez * r12
r22 = ez * r11
r23 = ct
!
! Loop on nodes (Voronoi centers) N0.
! LPL indexes the last neighbor of N0.
!
do n0 = 1, n
lpl = lend(n0)
!
! Set KV2 to the first (and last) vertex index and compute
! its coordinates (X2,Y2,Z2) in the rotated coordinate system.
!
kv2 = listc(lpl)
x2 = r11 * xc(kv2) + r12 * yc(kv2)
y2 = r21 * xc(kv2) + r22 * yc(kv2) + r23 * zc(kv2)
z2 = ex * xc(kv2) + ey * yc(kv2) + ez * zc(kv2)
!
! IN2 = TRUE iff KV2 is in the window.
!
in2 = ( 0.0D+00 <= z2 ) .and. ( x2 * x2 + y2 * y2 <= wrs )
!
! Loop on neighbors N1 of N0. For each triangulation edge
! N0-N1, KV1-KV2 is the corresponding Voronoi edge.
!
lp = lpl
do
lp = lptr(lp)
kv1 = kv2
x1 = x2
y1 = y2
z1 = z2
in1 = in2
kv2 = listc(lp)
!
! Compute the new values of (X2,Y2,Z2) and IN2.
!
x2 = r11 * xc(kv2) + r12 * yc(kv2)
y2 = r21 * xc(kv2) + r22 * yc(kv2) + r23 * zc(kv2)
z2 = ex * xc(kv2) + ey * yc(kv2) + ez * zc(kv2)
in2 = 0.0D+00 <= z2 .and. x2 * x2 + y2 * y2 <= wrs
!
! Add edge KV1-KV2 to the path iff both endpoints are inside
! the window and KV1 < KV2, or KV1 is inside and KV2 is
! outside (so that the edge is drawn only once).
!
if ( in1 .and. ( .not. in2 .or. kv1 < kv2 ) ) then
!
! If KV2 is a 'southern hemisphere' point, move it to the
! intersection of edge KV1-KV2 with the equator so that
! the edge is clipped properly. Z2 is implicitly set to 0.
!
if ( z2 < 0.0D+00 ) then
x2 = z1 * x2 - z2 * x1
y2 = z1 * y2 - z2 * y1
t = sqrt ( x2 * x2 + y2 * y2 )
x2 = x2 / t
y2 = y2 / t
end if
write ( lun, '(2f12.6,a,2f12.6,a)' ) &
x1, y1, ' moveto', x2, y2, ' lineto'
end if
if ( lp == lpl ) then
exit
end if
end do
end do
!
! Paint the path and restore the saved graphics state (with no clip path).
!
write ( lun, '(a)' ) 'stroke'
write ( lun, '(a)' ) 'grestore'
if ( numbr ) then
!
! Nodes in the window are to be labeled with their indexes.
! Convert FSIZN from points to world coordinates, and
! output the commands to select a font and scale it.
!
t = fsizn / sf
write ( lun, '(a)' ) '/Helvetica findfont'
write ( lun, '(f12.6,a)' ) t, ' scalefont setfont'
!
! Loop on visible nodes N0 that project to points (X0,Y0) in
! the window.
!
do n0 = 1, n
if ( ex * x(n0) + ey * y(n0) + ez * z(n0) < 0.0D+00 ) then
cycle
end if
x0 = r11 * x(n0) + r12 * y(n0)
y0 = r21 * x(n0) + r22 * y(n0) + r23 * z(n0)
!
! Move to (X0,Y0), and draw the label N0 with the origin
! of the first character at (X0,Y0).
!
if ( x0 * x0 + y0 * y0 <= wrs ) then
write ( lun, '(2f12.6,a)' ) x0, y0, ' moveto'
write ( lun, '(a,i3,a)' ) '(', n0, ') show'
end if
end do
end if
!
! Convert FSIZT from points to world coordinates, and output
! the commands to select a font and scale it.
!
t = fsizt / sf
write ( lun, '(a)' ) '/Helvetica findfont'
write ( lun, '(f12.6,a)' ) t, ' scalefont setfont'
!
! Display TITLE centered above the plot:
!
y0 = wr + 3.0D+00 * t
write ( lun, '(a)' ) title
write ( lun, '(a,g12.6,a)' ) ' stringwidth pop 2 div neg ', y0, ' moveto'
write ( lun, '(a)' ) title
write ( lun, '(a)' ) ' show'
!
! Display the window center and radius below the plot.
!
if ( annot ) then
x0 = -wr
y0 = -wr - 50.0D+00 / sf
write ( lun, '(2f12.6,a)' ) x0, y0, ' moveto'
write ( lun, '(a,f7.2,a,f8.2,a)' ) '(Window center: Latitude = ', elat, &
', Longitude = ', elon , ') show'
y0 = y0 - 2.0D+00 * t
write ( lun, '(2f12.6,a)' ) x0, y0, ' moveto'
write ( lun, '(a,f5.2,a)' ) '(Angular extent = ', a, ') show'
end if
!
! Paint the path and output the showpage command and end-of-file indicator.
!
write ( lun, '(a)' ) 'stroke'
write ( lun, '(a)' ) 'showpage'
write ( lun, '(a)' ) '%%eof'
return
end