SABER
tools_stripack.F90
Go to the documentation of this file.
1 !----------------------------------------------------------------------
2 ! Module: tools_stripack
3 !> STRIPACK routines
4 ! Source: https://people.sc.fsu.edu/~jburkardt/f_src/stripack/stripack.html
5 ! Author: Robert Renka
6 ! Original licensing: none
7 ! Modified by Benjamin Menetrier for BUMP
8 ! Licensing: this code is distributed under the CeCILL-C license
9 ! Copyright © 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
10 !----------------------------------------------------------------------
12 
13 use tools_kinds, only: kind_real
14 use tools_repro, only: eq,inf,sup,supeq,small
15 use type_mpl, only: mpl_type
16 
17 implicit none
18 
19 private
21 
22 contains
23 
24 subroutine addnod ( mpl, nst, k, x, y, z, list, lptr, lend, lnew, ier )
25 
26 !*****************************************************************************80
27 !
28 ! Subroutine: addnod
29 !> Add a node to a triangulation
30 !
31 ! Discussion:
32 !
33 ! This subroutine adds node K to a triangulation of the
34 ! convex hull of nodes 1, ..., K-1, producing a triangulation
35 ! of the convex hull of nodes 1, ..., K.
36 !
37 ! The algorithm consists of the following steps: node K
38 ! is located relative to the triangulation (TRFIND), its
39 ! index is added to the data structure (INTADD or BDYADD),
40 ! and a sequence of swaps (SWPTST and SWAP) are applied to
41 ! the arcs opposite K so that all arcs incident on node K
42 ! and opposite node K are locally optimal (satisfy the circumcircle test).
43 !
44 ! Thus, if a Delaunay triangulation of nodes 1 through K-1 is input,
45 ! a Delaunay triangulation of nodes 1 through K will be output.
46 !
47 ! Modified:
48 !
49 ! 15 May 2007
50 !
51 ! Author:
52 !
53 ! Robert Renka
54 !
55 ! Reference:
56 !
57 ! Robert Renka,
58 ! Algorithm 772: STRIPACK,
59 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
60 ! ACM Transactions on Mathematical Software,
61 ! Volume 23, Number 3, September 1997, pages 416-434.
62 !
63 ! Parameters:
64 !
65 ! Input, integer NST, the index of a node at which TRFIND
66 ! begins its search. Search time depends on the proximity of this node to
67 ! K. If NST < 1, the search is begun at node K-1.
68 !
69 ! Input, integer K, the nodal index (index for X, Y, Z, and
70 ! LEND) of the new node to be added. 4 <= K.
71 !
72 ! Input, real ( kind_real ) X(K), Y(K), Z(K), the coordinates of the nodes.
73 !
74 ! Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(K),
75 ! LNEW. On input, the data structure associated with the triangulation of
76 ! nodes 1 to K-1. On output, the data has been updated to include node
77 ! K. The array lengths are assumed to be large enough to add node K.
78 ! Refer to TRMESH.
79 !
80 ! Output, integer IER, error indicator:
81 ! 0 if no errors were encountered.
82 ! -1 if K is outside its valid range on input.
83 ! -2 if all nodes (including K) are collinear (lie on a common geodesic).
84 ! L if nodes L and K coincide for some L < K.
85 !
86 ! Local parameters:
87 !
88 ! B1,B2,B3 = Unnormalized barycentric coordinates returned by TRFIND.
89 ! I1,I2,I3 = Vertex indexes of a triangle containing K
90 ! IN1 = Vertex opposite K: first neighbor of IO2
91 ! that precedes IO1. IN1,IO1,IO2 are in
92 ! counterclockwise order.
93 ! IO1,IO2 = Adjacent neighbors of K defining an arc to
94 ! be tested for a swap
95 ! IST = Index of node at which TRFIND begins its search
96 ! KK = Local copy of K
97 ! KM1 = K-1
98 ! L = Vertex index (I1, I2, or I3) returned in IER
99 ! if node K coincides with a vertex
100 ! LP = LIST pointer
101 ! LPF = LIST pointer to the first neighbor of K
102 ! LPO1 = LIST pointer to IO1
103 ! LPO1S = Saved value of LPO1
104 ! P = Cartesian coordinates of node K
105 !
106  implicit none
107 
108  type(mpl_type),intent(inout) :: mpl !< MPI data
109 
110  integer k
111 
112  real ( kind_real ) b1
113  real ( kind_real ) b2
114  real ( kind_real ) b3
115  integer i1
116  integer i2
117  integer i3
118  integer ier
119  integer in1
120  integer io1
121  integer io2
122  integer ist
123  integer kk
124  integer km1
125  integer l
126  integer lend(k)
127  integer list(*)
128  integer lnew
129  integer lp
130  integer lpf
131  integer lpo1
132  integer lpo1s
133  integer lptr(*)
134  integer nst
135  real ( kind_real ) p(3)
136  logical output
137  real ( kind_real ) x(k)
138  real ( kind_real ) y(k)
139  real ( kind_real ) z(k)
140  character(len=1024),parameter :: subr = 'addnod'
141 
142  kk = k
143 
144  if ( kk < 4 ) then
145  ier = -1
146  write ( *, '(a)' ) ' '
147  write ( *, '(a)' ) 'ADDNOD - Fatal error!'
148  write ( *, '(a)' ) ' K < 4.'
149  call mpl%abort(subr,'stop in stripack')
150  end if
151 !
152 ! Initialization:
153 !
154  km1 = kk - 1
155  ist = nst
156  if ( ist < 1 ) then
157  ist = km1
158  end if
159 
160  p(1) = x(kk)
161  p(2) = y(kk)
162  p(3) = z(kk)
163 !
164 ! Find a triangle (I1,I2,I3) containing K or the rightmost
165 ! (I1) and leftmost (I2) visible boundary nodes as viewed
166 ! from node K.
167 !
168 
169  call trfind ( ist, p, km1, x, y, z, list, lptr, lend, b1, b2, b3, &
170  i1, i2, i3 )
171 !
172 ! Test for collinear or duplicate nodes.
173 !
174  if ( i1 == 0 ) call mpl%abort(subr,'colinear points')
175  if ( i3 /= 0 ) then
176 
177  l = i1
178 
179  if ( eq(p(1),x(l)) .and. eq(p(2),y(l)) .and. eq(p(3),z(l)) ) then
180  ier = l
181  return
182  end if
183 
184  l = i2
185 
186  if ( eq(p(1),x(l)) .and. eq(p(2),y(l)) .and. eq(p(3),z(l)) ) then
187  ier = l
188  return
189  end if
190 
191  l = i3
192  if ( eq(p(1),x(l)) .and. eq(p(2),y(l)) .and. eq(p(3),z(l)) ) then
193  ier = l
194  return
195  end if
196 
197  call intadd ( kk, i1, i2, i3, list, lptr, lend, lnew )
198 
199  else
200 
201  if ( i1 /= i2 ) then
202  call bdyadd ( kk, i1,i2, list, lptr, lend, lnew )
203  else
204  call covsph ( kk, i1, list, lptr, lend, lnew )
205  end if
206 
207  end if
208 
209  ier = 0
210 !
211 ! Initialize variables for optimization of the triangulation.
212 !
213  lp = lend(kk)
214  lpf = lptr(lp)
215  io2 = list(lpf)
216  lpo1 = lptr(lpf)
217  io1 = abs( list(lpo1) )
218 !
219 ! Begin loop: find the node opposite K.
220 !
221  do
222 
223  call lstptr ( lend(io1), io2, list, lptr, lp )
224 
225  if ( 0 <= list(lp) ) then
226 
227  lp = lptr(lp)
228  in1 = abs( list(lp) )
229 !
230 ! Swap test: if a swap occurs, two new arcs are
231 ! opposite K and must be tested.
232 !
233  lpo1s = lpo1
234 
235  call swptst ( in1, kk, io1, io2, x, y, z, output )
236  if ( .not. output ) then
237 
238  if ( lpo1 == lpf .or. list(lpo1) < 0 ) then
239  exit
240  end if
241 
242  io2 = io1
243  lpo1 = lptr(lpo1)
244  io1 = abs( list(lpo1) )
245  cycle
246 
247  end if
248 
249  call swap ( in1, kk, io1, io2, list, lptr, lend, lpo1 )
250 !
251 ! A swap is not possible because KK and IN1 are already
252 ! adjacent. This error in SWPTST only occurs in the
253 ! neutral case and when there are nearly duplicate nodes.
254 !
255  if ( lpo1 /= 0 ) then
256  io1 = in1
257  cycle
258  end if
259 
260  lpo1 = lpo1s
261 
262  end if
263 !
264 ! No swap occurred. Test for termination and reset IO2 and IO1.
265 !
266  if ( lpo1 == lpf .or. list(lpo1) < 0 ) then
267  exit
268  end if
269 
270  io2 = io1
271  lpo1 = lptr(lpo1)
272  io1 = abs( list(lpo1) )
273 
274  end do
275 
276  return
277 end subroutine addnod
278 subroutine bdyadd ( kk, i1, i2, list, lptr, lend, lnew )
279 
280 !*****************************************************************************80
281 !
282 ! Subroutine: bdyadd
283 !> Add a boundary node to a triangulation
284 !
285 ! Discussion:
286 !
287 ! This subroutine adds a boundary node to a triangulation
288 ! of a set of KK-1 points on the unit sphere. The data
289 ! structure is updated with the insertion of node KK, but no
290 ! optimization is performed.
291 !
292 ! This routine is identical to the similarly named routine
293 ! in TRIPACK.
294 !
295 ! Modified:
296 !
297 ! 16 June 2007
298 !
299 ! Author:
300 !
301 ! Robert Renka
302 !
303 ! Reference:
304 !
305 ! Robert Renka,
306 ! Algorithm 772: STRIPACK,
307 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
308 ! ACM Transactions on Mathematical Software,
309 ! Volume 23, Number 3, September 1997, pages 416-434.
310 !
311 ! Parameters:
312 !
313 ! Input, integer KK, the index of a node to be connected to
314 ! the sequence of all visible boundary nodes. 1 <= KK and
315 ! KK must not be equal to I1 or I2.
316 !
317 ! Input, integer I1, the first (rightmost as viewed from KK)
318 ! boundary node in the triangulation that is visible from
319 ! node KK (the line segment KK-I1 intersects no arcs.
320 !
321 ! Input, integer I2, the last (leftmost) boundary node that
322 ! is visible from node KK. I1 and I2 may be determined by TRFIND.
323 !
324 ! Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
325 ! LNEW, the triangulation data structure created by TRMESH.
326 ! Nodes I1 and I2 must be included
327 ! in the triangulation. On output, the data structure is updated with
328 ! the addition of node KK. Node KK is connected to I1, I2, and
329 ! all boundary nodes in between.
330 !
331 ! Local parameters:
332 !
333 ! K = Local copy of KK
334 ! LP = LIST pointer
335 ! LSAV = LIST pointer
336 ! N1,N2 = Local copies of I1 and I2, respectively
337 ! NEXT = Boundary node visible from K
338 ! NSAV = Boundary node visible from K
339 !
340  implicit none
341 
342  integer i1
343  integer i2
344  integer k
345  integer kk
346  integer lend(*)
347  integer list(*)
348  integer lnew
349  integer lp
350  integer lptr(*)
351  integer lsav
352  integer n1
353  integer n2
354  integer next
355  integer nsav
356 
357  k = kk
358  n1 = i1
359  n2 = i2
360 !
361 ! Add K as the last neighbor of N1.
362 !
363  lp = lend(n1)
364  lsav = lptr(lp)
365  lptr(lp) = lnew
366  list(lnew) = -k
367  lptr(lnew) = lsav
368  lend(n1) = lnew
369  lnew = lnew + 1
370  next = -list(lp)
371  list(lp) = next
372  nsav = next
373 !
374 ! Loop on the remaining boundary nodes between N1 and N2,
375 ! adding K as the first neighbor.
376 !
377  do
378 
379  lp = lend(next)
380  call insert ( k, lp, list, lptr, lnew )
381 
382  if ( next == n2 ) then
383  exit
384  end if
385 
386  next = -list(lp)
387  list(lp) = next
388 
389  end do
390 !
391 ! Add the boundary nodes between N1 and N2 as neighbors of node K.
392 !
393  lsav = lnew
394  list(lnew) = n1
395  lptr(lnew) = lnew + 1
396  lnew = lnew + 1
397  next = nsav
398 
399  do
400 
401  if ( next == n2 ) then
402  exit
403  end if
404 
405  list(lnew) = next
406  lptr(lnew) = lnew + 1
407  lnew = lnew + 1
408  lp = lend(next)
409  next = list(lp)
410 
411  end do
412 
413  list(lnew) = -n2
414  lptr(lnew) = lsav
415  lend(k) = lnew
416  lnew = lnew + 1
417 
418  return
419 end subroutine bdyadd
420 subroutine bnodes ( n, list, lptr, lend, nodes, nb, na, nt )
421 
422 !*****************************************************************************80
423 !
424 ! Subroutine: bnodes
425 !> Return the boundary nodes of a triangulation
426 !
427 ! Discussion:
428 !
429 ! Given a triangulation of N nodes on the unit sphere created by TRMESH,
430 ! this subroutine returns an array containing the indexes (if any) of
431 ! the counterclockwise sequence of boundary nodes, that is, the nodes on
432 ! the boundary of the convex hull of the set of nodes. The
433 ! boundary is empty if the nodes do not lie in a single
434 ! hemisphere. The numbers of boundary nodes, arcs, and
435 ! triangles are also returned.
436 !
437 ! Modified:
438 !
439 ! 16 June 2007
440 !
441 ! Author:
442 !
443 ! Robert Renka
444 !
445 ! Reference:
446 !
447 ! Robert Renka,
448 ! Algorithm 772: STRIPACK,
449 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
450 ! ACM Transactions on Mathematical Software,
451 ! Volume 23, Number 3, September 1997, pages 416-434.
452 !
453 ! Parameters:
454 !
455 ! Input, integer N, the number of nodes in the triangulation.
456 ! 3 <= N.
457 !
458 ! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), the
459 ! data structure defining the triangulation, created by TRMESH.
460 !
461 ! Output, integer NODES(*), the ordered sequence of NB boundary
462 ! node indexes in the range 1 to N. For safety, the dimension of NODES
463 ! should be N.
464 !
465 ! Output, integer NB, the number of boundary nodes.
466 !
467 ! Output, integer NA, NT, the number of arcs and triangles,
468 ! respectively, in the triangulation.
469 !
470 ! Local parameters:
471 !
472 ! K = NODES index
473 ! LP = LIST pointer
474 ! N0 = Boundary node to be added to NODES
475 ! NN = Local copy of N
476 ! NST = First element of nodes (arbitrarily chosen to be
477 ! the one with smallest index)
478 !
479  implicit none
480 
481  integer n
482 
483  integer i
484  integer k
485  integer lend(n)
486  integer list(6*(n-2))
487  integer lp
488  integer lptr(6*(n-2))
489  integer n0
490  integer na
491  integer nb
492  integer nn
493  integer nodes(*)
494  integer nst
495  integer nt
496 
497  nn = n
498 !
499 ! Search for a boundary node.
500 !
501  nst = 0
502 
503  do i = 1, nn
504 
505  lp = lend(i)
506 
507  if ( list(lp) < 0 ) then
508  nst = i
509  exit
510  end if
511 
512  end do
513 !
514 ! The triangulation contains no boundary nodes.
515 !
516  if ( nst == 0 ) then
517  nb = 0
518  na = 3 * ( nn - 2 )
519  nt = 2 * ( nn - 2 )
520  return
521  end if
522 !
523 ! NST is the first boundary node encountered.
524 !
525 ! Initialize for traversal of the boundary.
526 !
527  nodes(1) = nst
528  k = 1
529  n0 = nst
530 !
531 ! Traverse the boundary in counterclockwise order.
532 !
533  do
534 
535  lp = lend(n0)
536  lp = lptr(lp)
537  n0 = list(lp)
538 
539  if ( n0 == nst ) then
540  exit
541  end if
542 
543  k = k + 1
544  nodes(k) = n0
545 
546  end do
547 !
548 ! Store the counts.
549 !
550  nb = k
551  nt = 2 * n - nb - 2
552  na = nt + n - 1
553 
554  return
555 end subroutine bnodes
556 subroutine covsph ( kk, n0, list, lptr, lend, lnew )
557 
558 !*****************************************************************************80
559 !
560 ! Subroutine: covsph
561 !> Connect an exterior node to boundary nodes, covering the sphere
562 !
563 ! Discussion:
564 !
565 ! This subroutine connects an exterior node KK to all
566 ! boundary nodes of a triangulation of KK-1 points on the
567 ! unit sphere, producing a triangulation that covers the
568 ! sphere. The data structure is updated with the addition
569 ! of node KK, but no optimization is performed. All
570 ! boundary nodes must be visible from node KK.
571 !
572 ! Modified:
573 !
574 ! 16 June 2007
575 !
576 ! Author:
577 !
578 ! Robert Renka
579 !
580 ! Reference:
581 !
582 ! Robert Renka,
583 ! Algorithm 772: STRIPACK,
584 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
585 ! ACM Transactions on Mathematical Software,
586 ! Volume 23, Number 3, September 1997, pages 416-434.
587 !
588 ! Parameters:
589 !
590 ! Input, integer KK = Index of the node to be connected to the
591 ! set of all boundary nodes. 4 <= KK.
592 !
593 ! Input, integer N0 = Index of a boundary node (in the range
594 ! 1 to KK-1). N0 may be determined by TRFIND.
595 !
596 ! Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
597 ! LNEW, the triangulation data structure created by TRMESH. Node N0 must
598 ! be included in the triangulation. On output, updated with the addition
599 ! of node KK as the last entry. The updated triangulation contains no
600 ! boundary nodes.
601 !
602 ! Local parameters:
603 !
604 ! K = Local copy of KK
605 ! LP = LIST pointer
606 ! LSAV = LIST pointer
607 ! NEXT = Boundary node visible from K
608 ! NST = Local copy of N0
609 !
610  implicit none
611 
612  integer k
613  integer kk
614  integer lend(*)
615  integer list(*)
616  integer lnew
617  integer lp
618  integer lptr(*)
619  integer lsav
620  integer n0
621  integer next
622  integer nst
623 
624  k = kk
625  nst = n0
626 !
627 ! Traverse the boundary in clockwise order, inserting K as
628 ! the first neighbor of each boundary node, and converting
629 ! the boundary node to an interior node.
630 !
631  next = nst
632 
633  do
634 
635  lp = lend(next)
636  call insert ( k, lp, list, lptr, lnew )
637  next = -list(lp)
638  list(lp) = next
639 
640  if ( next == nst ) then
641  exit
642  end if
643 
644  end do
645 !
646 ! Traverse the boundary again, adding each node to K's adjacency list.
647 !
648  lsav = lnew
649 
650  do
651 
652  lp = lend(next)
653  list(lnew) = next
654  lptr(lnew) = lnew + 1
655  lnew = lnew + 1
656  next = list(lp)
657 
658  if ( next == nst ) then
659  exit
660  end if
661 
662  end do
663 
664  lptr(lnew-1) = lsav
665  lend(k) = lnew - 1
666 
667  return
668 end subroutine covsph
669 subroutine det ( x1, y1, z1, x2, y2, z2, x0, y0, z0, output )
670 !*****************************************************************************80
671 !
672 ! Subroutine: det
673 !> Compute 3D determinant
674 !
675 ! DET(X1,...,Z0) >= 0 if and only if (X0,Y0,Z0) is in the
676 ! (closed) left hemisphere defined by the plane containing (0,0,0),
677 ! (X1,Y1,Z1), and (X2,Y2,Z2), where left is defined relative to an
678 ! observer at (X1,Y1,Z1) facing (X2,Y2,Z2).
679 !
680  implicit none
681 
682  real ( kind_real ) x1
683  real ( kind_real ) y1
684  real ( kind_real ) z1
685  real ( kind_real ) x2
686  real ( kind_real ) y2
687  real ( kind_real ) z2
688  real ( kind_real ) x0
689  real ( kind_real ) y0
690  real ( kind_real ) z0
691  real ( kind_real ) t1
692  real ( kind_real ) t2
693  real ( kind_real ) t3
694  real ( kind_real ) output
695 
696  t1 = x0*(y1*z2-y2*z1)
697  t2 = y0*(x1*z2-x2*z1)
698  t3 = z0*(x1*y2-x2*y1)
699  output = t1 - t2 + t3
700 
701  ! Indistinguishability threshold for cross-plateform reproducibility
702  if (small(output,t1).or.small(output,t2).or.small(output,t3)) output = 0.0
703 end subroutine det
704 subroutine insert ( k, lp, list, lptr, lnew )
705 
706 !*****************************************************************************80
707 !
708 ! Subroutine: insert
709 !> Insert K as a neighbor of N1
710 !
711 ! Discussion:
712 !
713 ! This subroutine inserts K as a neighbor of N1 following
714 ! N2, where LP is the LIST pointer of N2 as a neighbor of
715 ! N1. Note that, if N2 is the last neighbor of N1, K will
716 ! become the first neighbor (even if N1 is a boundary node).
717 !
718 ! This routine is identical to the similarly named routine in TRIPACK.
719 !
720 ! Modified:
721 !
722 ! 16 June 2007
723 !
724 ! Author:
725 !
726 ! Robert Renka
727 !
728 ! Reference:
729 !
730 ! Robert Renka,
731 ! Algorithm 772: STRIPACK,
732 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
733 ! ACM Transactions on Mathematical Software,
734 ! Volume 23, Number 3, September 1997, pages 416-434.
735 !
736 ! Parameters:
737 !
738 ! Input, integer K, the index of the node to be inserted.
739 !
740 ! Input, integer LP, the LIST pointer of N2 as a neighbor of N1.
741 !
742 ! Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LNEW,
743 ! the data structure defining the triangulation, created by TRMESH.
744 ! On output, updated with the addition of node K.
745 !
746  implicit none
747 
748  integer k
749  integer list(*)
750  integer lnew
751  integer lp
752  integer lptr(*)
753  integer lsav
754 
755  lsav = lptr(lp)
756  lptr(lp) = lnew
757  list(lnew) = k
758  lptr(lnew) = lsav
759  lnew = lnew + 1
760 
761  return
762 end subroutine insert
763 function inside ( p, lv, xv, yv, zv, nv, listv, ier )
764 
765 !*****************************************************************************80
766 !
767 ! Function: inside
768 !> Determine if a point is inside a polygonal region
769 !
770 ! Discussion:
771 !
772 ! This function locates a point P relative to a polygonal
773 ! region R on the surface of the unit sphere, returning
774 ! INSIDE = TRUE if and only if P is contained in R. R is
775 ! defined by a cyclically ordered sequence of vertices which
776 ! form a positively-oriented simple closed curve. Adjacent
777 ! vertices need not be distinct but the curve must not be
778 ! self-intersecting. Also, while polygon edges are by definition
779 ! restricted to a single hemisphere, R is not so
780 ! restricted. Its interior is the region to the left as the
781 ! vertices are traversed in order.
782 !
783 ! The algorithm consists of selecting a point Q in R and
784 ! then finding all points at which the great circle defined
785 ! by P and Q intersects the boundary of R. P lies inside R
786 ! if and only if there is an even number of intersection
787 ! points between Q and P. Q is taken to be a point immediately
788 ! to the left of a directed boundary edge -- the first
789 ! one that results in no consistency-check failures.
790 !
791 ! If P is close to the polygon boundary, the problem is
792 ! ill-conditioned and the decision may be incorrect. Also,
793 ! an incorrect decision may result from a poor choice of Q
794 ! (if, for example, a boundary edge lies on the great circle
795 ! defined by P and Q). A more reliable result could be
796 ! obtained by a sequence of calls to INSIDE with the vertices
797 ! cyclically permuted before each call (to alter the
798 ! choice of Q).
799 !
800 ! Modified:
801 !
802 ! 16 June 2007
803 !
804 ! Author:
805 !
806 ! Robert Renka
807 !
808 ! Reference:
809 !
810 ! Robert Renka,
811 ! Algorithm 772: STRIPACK,
812 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
813 ! ACM Transactions on Mathematical Software,
814 ! Volume 23, Number 3, September 1997, pages 416-434.
815 !
816 ! Parameters:
817 !
818 ! Input, real ( kind_real ) P(3), the coordinates of the point (unit vector)
819 ! to be located.
820 !
821 ! Input, integer LV, the length of arrays XV, YV, and ZV.
822 !
823 ! Input, real ( kind_real ) XV(LV), YV(LV), ZV(LV), the coordinates of unit
824 ! vectors (points on the unit sphere).
825 !
826 ! Input, integer NV, the number of vertices in the polygon.
827 ! 3 <= NV <= LV.
828 !
829 ! Input, integer LISTV(NV), the indexes (for XV, YV, and ZV)
830 ! of a cyclically-ordered (and CCW-ordered) sequence of vertices that
831 ! define R. The last vertex (indexed by LISTV(NV)) is followed by the
832 ! first (indexed by LISTV(1)). LISTV entries must be in the range 1 to LV.
833 !
834 ! Output, logical INSIDE, TRUE if and only if P lies inside R unless
835 ! IER /= 0, in which case the value is not altered.
836 !
837 ! Output, integer IER, error indicator:
838 ! 0, if no errors were encountered.
839 ! 1, if LV or NV is outside its valid range.
840 ! 2, if a LISTV entry is outside its valid range.
841 ! 3, if the polygon boundary was found to be self-intersecting. This
842 ! error will not necessarily be detected.
843 ! 4, if every choice of Q (one for each boundary edge) led to failure of
844 ! some internal consistency check. The most likely cause of this error
845 ! is invalid input: P = (0,0,0), a null or self-intersecting polygon, etc.
846 !
847 ! Local parameters:
848 !
849 ! B = Intersection point between the boundary and
850 ! the great circle defined by P and Q.
851 !
852 ! BP,BQ = <B,P> and <B,Q>, respectively, maximized over
853 ! intersection points B that lie between P and
854 ! Q (on the shorter arc) -- used to find the
855 ! closest intersection points to P and Q
856 ! CN = Q X P = normal to the plane of P and Q
857 ! D = Dot product <B,P> or <B,Q>
858 ! EPS = Parameter used to define Q as the point whose
859 ! orthogonal distance to (the midpoint of)
860 ! boundary edge V1->V2 is approximately EPS/
861 ! (2*Cos(A/2)), where <V1,V2> = Cos(A).
862 ! EVEN = TRUE iff an even number of intersection points
863 ! lie between P and Q (on the shorter arc)
864 ! I1,I2 = Indexes (LISTV elements) of a pair of adjacent
865 ! boundary vertices (endpoints of a boundary
866 ! edge)
867 ! IERR = Error flag for calls to INTRSC (not tested)
868 ! IMX = Local copy of LV and maximum value of I1 and I2
869 ! K = DO-loop index and LISTV index
870 ! K0 = LISTV index of the first endpoint of the
871 ! boundary edge used to compute Q
872 ! LFT1,LFT2 = Logical variables associated with I1 and I2 in
873 ! the boundary traversal: TRUE iff the vertex
874 ! is strictly to the left of Q->P (<V,CN> > 0)
875 ! N = Local copy of NV
876 ! NI = Number of intersections (between the boundary
877 ! curve and the great circle P-Q) encountered
878 ! PINR = TRUE iff P is to the left of the directed
879 ! boundary edge associated with the closest
880 ! intersection point to P that lies between P
881 ! and Q (a left-to-right intersection as
882 ! viewed from Q), or there is no intersection
883 ! between P and Q (on the shorter arc)
884 ! PN,QN = P X CN and CN X Q, respectively: used to
885 ! locate intersections B relative to arc Q->P
886 ! Q = (V1 + V2 + EPS*VN/VNRM)/QNRM, where V1->V2 is
887 ! the boundary edge indexed by LISTV(K0) ->
888 ! LISTV(K0+1)
889 ! QINR = TRUE iff Q is to the left of the directed
890 ! boundary edge associated with the closest
891 ! intersection point to Q that lies between P
892 ! and Q (a right-to-left intersection as
893 ! viewed from Q), or there is no intersection
894 ! between P and Q (on the shorter arc)
895 ! QNRM = Euclidean norm of V1+V2+EPS*VN/VNRM used to
896 ! compute (normalize) Q
897 ! V1,V2 = Vertices indexed by I1 and I2 in the boundary
898 ! traversal
899 ! VN = V1 X V2, where V1->V2 is the boundary edge
900 ! indexed by LISTV(K0) -> LISTV(K0+1)
901 ! VNRM = Euclidean norm of VN
902 !
903  implicit none
904 
905  integer lv
906  integer nv
907 
908  real ( kind_real ) b(3)
909  real ( kind_real ) bp
910  real ( kind_real ) bq
911  real ( kind_real ) cn(3)
912  real ( kind_real ) d
913  real ( kind_real ), parameter :: eps = 0.001_kind_real
914  logical even
915  integer i1
916  integer i2
917  integer ier
918  integer ierr
919  integer imx
920  logical inside
921  integer k
922  integer k0
923  logical lft1
924  logical lft2
925  integer listv(nv)
926  integer n
927  integer ni
928  real ( kind_real ) p(3)
929  logical pinr
930  real ( kind_real ) pn(3)
931  real ( kind_real ) q(3)
932  logical qinr
933  real ( kind_real ) qn(3)
934  real ( kind_real ) qnrm
935  real ( kind_real ) v1(3)
936  real ( kind_real ) v2(3)
937  real ( kind_real ) vn(3)
938  real ( kind_real ) vnrm
939  real ( kind_real ) xv(lv)
940  real ( kind_real ) yv(lv)
941  real ( kind_real ) zv(lv)
942 !
943 ! Default value
944 !
945  inside = .false.
946 !
947 ! Store local parameters.
948 !
949  imx = lv
950  n = nv
951 !
952 ! Test for error 1.
953 !
954  if ( n < 3 .or. imx < n ) then
955  ier = 1
956  return
957  end if
958 !
959 ! Initialize K0.
960 !
961  k0 = 0
962  i1 = listv(1)
963 
964  if ( i1 < 1 .or. imx < i1 ) then
965  ier = 2
966  return
967  end if
968 !
969 ! Increment K0 and set Q to a point immediately to the left
970 ! of the midpoint of edge V1->V2 = LISTV(K0)->LISTV(K0+1):
971 ! Q = (V1 + V2 + EPS*VN/VNRM)/QNRM, where VN = V1 X V2.
972 !
973 1 continue
974 
975  k0 = k0 + 1
976 
977  if ( n < k0 ) then
978  ier = 4
979  return
980  end if
981 
982  i1 = listv(k0)
983 
984  if ( k0 < n ) then
985  i2 = listv(k0+1)
986  else
987  i2 = listv(1)
988  end if
989 
990  if ( i2 < 1 .or. imx < i2 ) then
991  ier = 2
992  return
993  end if
994 
995  vn(1) = yv(i1) * zv(i2) - zv(i1) * yv(i2)
996  vn(2) = zv(i1) * xv(i2) - xv(i1) * zv(i2)
997  vn(3) = xv(i1) * yv(i2) - yv(i1) * xv(i2)
998  vnrm = sqrt( sum( vn(1:3)**2 ) )
999 
1000  if ( eq(vnrm,0.0_kind_real) ) then
1001  go to 1
1002  end if
1003 
1004  q(1) = xv(i1) + xv(i2) + eps * vn(1) / vnrm
1005  q(2) = yv(i1) + yv(i2) + eps * vn(2) / vnrm
1006  q(3) = zv(i1) + zv(i2) + eps * vn(3) / vnrm
1007 
1008  qnrm = sqrt( sum( q(1:3)**2 ) )
1009 
1010  q(1) = q(1) / qnrm
1011  q(2) = q(2) / qnrm
1012  q(3) = q(3) / qnrm
1013 !
1014 ! Compute CN = Q X P, PN = P X CN, and QN = CN X Q.
1015 !
1016  cn(1) = q(2) * p(3) - q(3) * p(2)
1017  cn(2) = q(3) * p(1) - q(1) * p(3)
1018  cn(3) = q(1) * p(2) - q(2) * p(1)
1019 
1020  if ( eq(cn(1),0.0_kind_real) .and. eq(cn(2),0.0_kind_real) .and. eq(cn(2),0.0_kind_real) ) then
1021  go to 1
1022  end if
1023 
1024  pn(1) = p(2) * cn(3) - p(3) * cn(2)
1025  pn(2) = p(3) * cn(1) - p(1) * cn(3)
1026  pn(3) = p(1) * cn(2) - p(2) * cn(1)
1027  qn(1) = cn(2) * q(3) - cn(3) * q(2)
1028  qn(2) = cn(3) * q(1) - cn(1) * q(3)
1029  qn(3) = cn(1) * q(2) - cn(2) * q(1)
1030 !
1031 ! Initialize parameters for the boundary traversal.
1032 !
1033  ni = 0
1034  even = .true.
1035  bp = -2.0_kind_real
1036  bq = -2.0_kind_real
1037  pinr = .true.
1038  qinr = .true.
1039  i2 = listv(n)
1040 
1041  if ( i2 < 1 .or. imx < i2 ) then
1042  ier = 2
1043  return
1044  end if
1045 
1046  lft2 = sup(cn(1) * xv(i2) + cn(2) * yv(i2) + cn(3) * zv(i2),0.0_kind_real)
1047 !
1048 ! Loop on boundary arcs I1->I2.
1049 !
1050  do k = 1, n
1051 
1052  i1 = i2
1053  lft1 = lft2
1054  i2 = listv(k)
1055 
1056  if ( i2 < 1 .or. imx < i2 ) then
1057  ier = 2
1058  return
1059  end if
1060 
1061  lft2 = sup(cn(1) * xv(i2) + cn(2) * yv(i2) + cn(3) * zv(i2),0.0_kind_real)
1062 
1063  if ( lft1 .eqv. lft2 ) then
1064  cycle
1065  end if
1066 !
1067 ! I1 and I2 are on opposite sides of Q->P. Compute the
1068 ! point of intersection B.
1069 !
1070  ni = ni + 1
1071  v1(1) = xv(i1)
1072  v1(2) = yv(i1)
1073  v1(3) = zv(i1)
1074  v2(1) = xv(i2)
1075  v2(2) = yv(i2)
1076  v2(3) = zv(i2)
1077 
1078  call intrsc ( v1, v2, cn, b, ierr )
1079 !
1080 ! B is between Q and P (on the shorter arc) iff
1081 ! B Forward Q->P and B Forward P->Q iff
1082 ! <B,QN> > 0 and 0 < <B,PN>.
1083 !
1084  if ( sup(dot_product( b(1:3), qn(1:3) ),0.0_kind_real) .and. &
1085  sup(dot_product( b(1:3), pn(1:3) ),0.0_kind_real) ) then
1086 !
1087 ! Update EVEN, BQ, QINR, BP, and PINR.
1088 !
1089  even = .not. even
1090  d = dot_product( b(1:3), q(1:3) )
1091 
1092  if ( sup(d,bq) ) then
1093  bq = d
1094  qinr = lft2
1095  end if
1096 
1097  d = dot_product( b(1:3), p(1:3) )
1098 
1099  if ( sup(d,bq) ) then
1100  bp = d
1101  pinr = lft1
1102  end if
1103 
1104  end if
1105 
1106  end do
1107 !
1108 ! Test for consistency: NI must be even and QINR must be TRUE.
1109 !
1110  if ( ni /= 2 * ( ni / 2 ) .or. .not. qinr ) then
1111  go to 1
1112  end if
1113 !
1114 ! Test for error 3: different values of PINR and EVEN.
1115 !
1116  if ( pinr .neqv. even ) then
1117  ier = 3
1118  return
1119  end if
1120 
1121  ier = 0
1122  inside = even
1123 
1124  return
1125 end function inside
1126 subroutine intadd ( kk, i1, i2, i3, list, lptr, lend, lnew )
1127 
1128 !*****************************************************************************80
1129 !
1130 ! Subroutine: intadd
1131 !> Add an interior node to a triangulation
1132 !
1133 ! Discussion:
1134 !
1135 ! This subroutine adds an interior node to a triangulation
1136 ! of a set of points on the unit sphere. The data structure
1137 ! is updated with the insertion of node KK into the triangle
1138 ! whose vertices are I1, I2, and I3. No optimization of the
1139 ! triangulation is performed.
1140 !
1141 ! This routine is identical to the similarly named routine in TRIPACK.
1142 !
1143 ! Modified:
1144 !
1145 ! 16 June 2007
1146 !
1147 ! Author:
1148 !
1149 ! Robert Renka
1150 !
1151 ! Reference:
1152 !
1153 ! Robert Renka,
1154 ! Algorithm 772: STRIPACK,
1155 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
1156 ! ACM Transactions on Mathematical Software,
1157 ! Volume 23, Number 3, September 1997, pages 416-434.
1158 !
1159 ! Parameters:
1160 !
1161 ! Input, integer KK, the index of the node to be inserted.
1162 ! 1 <= KK and KK must not be equal to I1, I2, or I3.
1163 !
1164 ! Input, integer I1, I2, I3, indexes of the
1165 ! counterclockwise-ordered sequence of vertices of a triangle which contains
1166 ! node KK.
1167 !
1168 ! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), LNEW,
1169 ! the data structure defining the triangulation, created by TRMESH. Triangle
1170 ! (I1,I2,I3) must be included in the triangulation.
1171 ! On output, updated with the addition of node KK. KK
1172 ! will be connected to nodes I1, I2, and I3.
1173 !
1174 ! Local parameters:
1175 !
1176 ! K = Local copy of KK
1177 ! LP = LIST pointer
1178 ! N1,N2,N3 = Local copies of I1, I2, and I3
1179 !
1180  implicit none
1181 
1182  integer i1
1183  integer i2
1184  integer i3
1185  integer k
1186  integer kk
1187  integer lend(*)
1188  integer list(*)
1189  integer lnew
1190  integer lp
1191  integer lptr(*)
1192  integer n1
1193  integer n2
1194  integer n3
1195 
1196  k = kk
1197 !
1198 ! Initialization.
1199 !
1200  n1 = i1
1201  n2 = i2
1202  n3 = i3
1203 !
1204 ! Add K as a neighbor of I1, I2, and I3.
1205 !
1206  call lstptr ( lend(n1), n2, list, lptr, lp )
1207  call insert ( k, lp, list, lptr, lnew )
1208 
1209  call lstptr ( lend(n2), n3, list, lptr, lp )
1210  call insert ( k, lp, list, lptr, lnew )
1211 
1212  call lstptr ( lend(n3), n1, list, lptr, lp )
1213  call insert ( k, lp, list, lptr, lnew )
1214 !
1215 ! Add I1, I2, and I3 as neighbors of K.
1216 !
1217  list(lnew) = n1
1218  list(lnew+1) = n2
1219  list(lnew+2) = n3
1220  lptr(lnew) = lnew + 1
1221  lptr(lnew+1) = lnew + 2
1222  lptr(lnew+2) = lnew
1223  lend(k) = lnew + 2
1224  lnew = lnew + 3
1225 
1226  return
1227 end subroutine intadd
1228 subroutine intrsc ( p1, p2, cn, p, ier )
1229 
1230 !*****************************************************************************80
1231 !
1232 ! Subroutine: intrsc
1233 !> Find the intersection of two great circles
1234 !
1235 ! Discussion:
1236 !
1237 ! Given a great circle C and points P1 and P2 defining an
1238 ! arc A on the surface of the unit sphere, where A is the
1239 ! shorter of the two portions of the great circle C12
1240 ! associated with P1 and P2, this subroutine returns the point
1241 ! of intersection P between C and C12 that is closer to A.
1242 ! Thus, if P1 and P2 lie in opposite hemispheres defined by
1243 ! C, P is the point of intersection of C with A.
1244 !
1245 ! Modified:
1246 !
1247 ! 16 June 2007
1248 !
1249 ! Author:
1250 !
1251 ! Robert Renka
1252 !
1253 ! Reference:
1254 !
1255 ! Robert Renka,
1256 ! Algorithm 772: STRIPACK,
1257 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
1258 ! ACM Transactions on Mathematical Software,
1259 ! Volume 23, Number 3, September 1997, pages 416-434.
1260 !
1261 ! Parameters:
1262 !
1263 ! Input, real ( kind_real ) P1(3), P2(3), the coordinates of unit vectors.
1264 !
1265 ! Input, real ( kind_real ) CN(3), the coordinates of a nonzero vector
1266 ! which defines C as the intersection of the plane whose normal is CN
1267 ! with the unit sphere. Thus, if C is to be the great circle defined
1268 ! by P and Q, CN should be P X Q.
1269 !
1270 ! Output, real ( kind_real ) P(3), point of intersection defined above
1271 ! unless IER is not 0, in which case P is not altered.
1272 !
1273 ! Output, integer IER, error indicator.
1274 ! 0, if no errors were encountered.
1275 ! 1, if <CN,P1> = <CN,P2>. This occurs iff P1 = P2 or CN = 0 or there are
1276 ! two intersection points at the same distance from A.
1277 ! 2, if P2 = -P1 and the definition of A is therefore ambiguous.
1278 !
1279 ! Local parameters:
1280 !
1281 ! D1 = <CN,P1>
1282 ! D2 = <CN,P2>
1283 ! I = DO-loop index
1284 ! PP = P1 + T*(P2-P1) = Parametric representation of the
1285 ! line defined by P1 and P2
1286 ! PPN = Norm of PP
1287 ! T = D1/(D1-D2) = Parameter value chosen so that PP lies
1288 ! in the plane of C
1289 !
1290  implicit none
1291 
1292  real ( kind_real ) cn(3)
1293  real ( kind_real ) d1
1294  real ( kind_real ) d2
1295  integer ier
1296  real ( kind_real ) p(3)
1297  real ( kind_real ) p1(3)
1298  real ( kind_real ) p2(3)
1299  real ( kind_real ) pp(3)
1300  real ( kind_real ) ppn
1301  real ( kind_real ) t
1302 
1303  d1 = dot_product( cn(1:3), p1(1:3) )
1304  d2 = dot_product( cn(1:3), p2(1:3) )
1305 
1306  if ( eq(d1,d2) ) then
1307  ier = 1
1308  return
1309  end if
1310 !
1311 ! Solve for T such that <PP,CN> = 0 and compute PP and PPN.
1312 !
1313  t = d1 / ( d1 - d2 )
1314 
1315  pp(1:3) = p1(1:3) + t * ( p2(1:3) - p1(1:3) )
1316 
1317  ppn = dot_product( pp(1:3), pp(1:3) )
1318 !
1319 ! PPN = 0 iff PP = 0 iff P2 = -P1 (and T = .5).
1320 !
1321  if ( eq(ppn,0.0_kind_real) ) then
1322  ier = 2
1323  return
1324  end if
1325 
1326  ppn = sqrt( ppn )
1327 !
1328 ! Compute P = PP/PPN.
1329 !
1330  p(1:3) = pp(1:3) / ppn
1331 
1332  ier = 0
1333 
1334  return
1335 end subroutine intrsc
1336 subroutine jrand ( n, ix, iy, iz, output )
1337 
1338 !*****************************************************************************80
1339 !
1340 ! Subroutine: jrand
1341 !> Return a random integer between 1 and N
1342 !
1343 ! Discussion:
1344 !
1345 ! This function returns a uniformly distributed pseudorandom integer
1346 ! in the range 1 to N.
1347 !
1348 ! Modified:
1349 !
1350 ! 16 June 2007
1351 !
1352 ! Author:
1353 !
1354 ! Robert Renka
1355 !
1356 ! Reference:
1357 !
1358 ! Brian Wichmann, David Hill,
1359 ! An Efficient and Portable Pseudo-random Number Generator,
1360 ! Applied Statistics,
1361 ! Volume 31, Number 2, 1982, pages 188-190.
1362 !
1363 ! Parameters:
1364 !
1365 ! Input, integer N, the maximum value to be returned.
1366 !
1367 ! Input/output, integer IX, IY, IZ = seeds initialized to
1368 ! values in the range 1 to 30,000 before the first call to JRAND, and
1369 ! not altered between subsequent calls (unless a sequence of random
1370 ! numbers is to be repeated by reinitializing the seeds).
1371 !
1372 ! Output, integer JRAND, a random integer in the range 1 to N.
1373 !
1374 ! Local parameters:
1375 !
1376 ! U = Pseudo-random number uniformly distributed in the interval (0,1).
1377 ! X = Pseudo-random number in the range 0 to 3 whose fractional part is U.
1378 !
1379  implicit none
1380 
1381  integer ix
1382  integer iy
1383  integer iz
1384  integer output
1385  integer n
1386  real ( kind_real ) u
1387  real ( kind_real ) x
1388 
1389  ix = mod( 171 * ix, 30269 )
1390  iy = mod( 172 * iy, 30307 )
1391  iz = mod( 170 * iz, 30323 )
1392 
1393  x = ( real( ix, kind_real ) / 30269.0_kind_real ) &
1394  + ( real( iy, kind_real ) / 30307.0_kind_real ) &
1395  + ( real( iz, kind_real ) / 30323.0_kind_real )
1396 
1397  u = x - int( x )
1398  output = int( real( n, kind_real ) * u ) + 1
1399 
1400  return
1401 end subroutine jrand
1402 subroutine left ( x1, y1, z1, x2, y2, z2, x0, y0, z0, output )
1403 
1404 !*****************************************************************************80
1405 !
1406 ! Subroutine: left
1407 !> Determin whether a node is to the left of a plane through the origin
1408 !
1409 ! Discussion:
1410 !
1411 ! This function determines whether node N0 is in the
1412 ! (closed) left hemisphere defined by the plane containing
1413 ! N1, N2, and the origin, where left is defined relative to
1414 ! an observer at N1 facing N2.
1415 !
1416 ! Modified:
1417 !
1418 ! 16 June 2007
1419 !
1420 ! Author:
1421 !
1422 ! Robert Renka
1423 !
1424 ! Reference:
1425 !
1426 ! Robert Renka,
1427 ! Algorithm 772: STRIPACK,
1428 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
1429 ! ACM Transactions on Mathematical Software,
1430 ! Volume 23, Number 3, September 1997, pages 416-434.
1431 !
1432 ! Parameters:
1433 !
1434 ! Input, real ( kind_real ) X1, Y1, Z1 = Coordinates of N1.
1435 !
1436 ! Input, real ( kind_real ) X2, Y2, Z2 = Coordinates of N2.
1437 !
1438 ! Input, real ( kind_real ) X0, Y0, Z0 = Coordinates of N0.
1439 !
1440 ! Output, logical LEFT = TRUE if and only if N0 is in the closed
1441 ! left hemisphere.
1442 !
1443  implicit none
1444 
1445  logical output
1446  real ( kind_real ) x0
1447  real ( kind_real ) x1
1448  real ( kind_real ) x2
1449  real ( kind_real ) y0
1450  real ( kind_real ) y1
1451  real ( kind_real ) y2
1452  real ( kind_real ) z0
1453  real ( kind_real ) z1
1454  real ( kind_real ) z2
1455  real ( kind_real ) zz
1456 !
1457 ! LEFT = TRUE iff <N0,N1 X N2> = det(N0,N1,N2) >= 0.
1458 !
1459 
1460  call det ( x1, y1, z1, x2, y2, z2, x0, y0, z0, zz )
1461  output = sup(zz,0.0_kind_real)
1462 
1463  return
1464 end subroutine left
1465 subroutine lstptr ( lpl, nb, list, lptr, output )
1466 
1467 !*****************************************************************************80
1468 !
1469 ! Subroutine: lstptr
1470 !> Return the index of NB in the adjacency list
1471 !
1472 ! Discussion:
1473 !
1474 ! This function returns the index (LIST pointer) of NB in
1475 ! the adjacency list for N0, where LPL = LEND(N0).
1476 !
1477 ! This function is identical to the similarly named function in TRIPACK.
1478 !
1479 ! Modified:
1480 !
1481 ! 16 June 2007
1482 !
1483 ! Author:
1484 !
1485 ! Robert Renka
1486 !
1487 ! Reference:
1488 !
1489 ! Robert Renka,
1490 ! Algorithm 772: STRIPACK,
1491 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
1492 ! ACM Transactions on Mathematical Software,
1493 ! Volume 23, Number 3, September 1997, pages 416-434.
1494 !
1495 ! Parameters:
1496 !
1497 ! Input, integer LPL, is LEND(N0).
1498 !
1499 ! Input, integer NB, index of the node whose pointer is to
1500 ! be returned. NB must be connected to N0.
1501 !
1502 ! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), the data
1503 ! structure defining the triangulation, created by TRMESH.
1504 !
1505 ! Output, integer LSTPTR, pointer such that LIST(LSTPTR) = NB or
1506 ! LIST(LSTPTR) = -NB, unless NB is not a neighbor of N0, in which
1507 ! case LSTPTR = LPL.
1508 !
1509 ! Local parameters:
1510 !
1511 ! LP = LIST pointer
1512 ! ND = Nodal index
1513 !
1514  implicit none
1515 
1516  integer list(*)
1517  integer lp
1518  integer lpl
1519  integer lptr(*)
1520  integer output
1521  integer nb
1522  integer nd
1523 
1524  lp = lptr(lpl)
1525 
1526  do
1527 
1528  nd = list(lp)
1529 
1530  if ( nd == nb ) then
1531  exit
1532  end if
1533 
1534  lp = lptr(lp)
1535 
1536  if ( lp == lpl ) then
1537  exit
1538  end if
1539 
1540  end do
1541 
1542  output = lp
1543 
1544  return
1545 end subroutine lstptr
1546 subroutine swap ( in1, in2, io1, io2, list, lptr, lend, lp21 )
1547 
1548 !*****************************************************************************80
1549 !
1550 ! Subroutine: swap
1551 !> Replace the diagonal arc of a quadrilateral with the other diagonal
1552 !
1553 ! Discussion:
1554 !
1555 ! Given a triangulation of a set of points on the unit
1556 ! sphere, this subroutine replaces a diagonal arc in a
1557 ! strictly convex quadrilateral (defined by a pair of adja-
1558 ! cent triangles) with the other diagonal. Equivalently, a
1559 ! pair of adjacent triangles is replaced by another pair
1560 ! having the same union.
1561 !
1562 ! Modified:
1563 !
1564 ! 16 June 2007
1565 !
1566 ! Author:
1567 !
1568 ! Robert Renka
1569 !
1570 ! Reference:
1571 !
1572 ! Robert Renka,
1573 ! Algorithm 772: STRIPACK,
1574 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
1575 ! ACM Transactions on Mathematical Software,
1576 ! Volume 23, Number 3, September 1997, pages 416-434.
1577 !
1578 ! Parameters:
1579 !
1580 ! Input, integer IN1, IN2, IO1, IO2, nodal indexes of the
1581 ! vertices of the quadrilateral. IO1-IO2 is replaced by IN1-IN2.
1582 ! (IO1,IO2,IN1) and (IO2,IO1,IN2) must be triangles on input.
1583 !
1584 ! Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
1585 ! the data structure defining the triangulation, created by TRMESH.
1586 ! On output, updated with the swap; triangles (IO1,IO2,IN1) an (IO2,IO1,IN2)
1587 ! are replaced by (IN1,IN2,IO2) and (IN2,IN1,IO1) unless LP21 = 0.
1588 !
1589 ! Output, integer LP21, index of IN1 as a neighbor of IN2 after
1590 ! the swap is performed unless IN1 and IN2 are adjacent on input, in which
1591 ! case LP21 = 0.
1592 !
1593 ! Local parameters:
1594 !
1595 ! LP, LPH, LPSAV = LIST pointers
1596 !
1597  implicit none
1598 
1599  integer in1
1600  integer in2
1601  integer io1
1602  integer io2
1603  integer lend(*)
1604  integer list(*)
1605  integer lp
1606  integer lp21
1607  integer lph
1608  integer lpsav
1609  integer lptr(*)
1610 !
1611 ! Test for IN1 and IN2 adjacent.
1612 !
1613  call lstptr ( lend(in1), in2, list, lptr, lp )
1614 
1615  if ( abs( list(lp) ) == in2 ) then
1616  lp21 = 0
1617  return
1618  end if
1619 !
1620 ! Delete IO2 as a neighbor of IO1.
1621 !
1622  call lstptr ( lend(io1), in2, list, lptr, lp )
1623  lph = lptr(lp)
1624  lptr(lp) = lptr(lph)
1625 !
1626 ! If IO2 is the last neighbor of IO1, make IN2 the last neighbor.
1627 !
1628  if ( lend(io1) == lph ) then
1629  lend(io1) = lp
1630  end if
1631 !
1632 ! Insert IN2 as a neighbor of IN1 following IO1 using the hole created above.
1633 !
1634  call lstptr ( lend(in1), io1, list, lptr, lp )
1635  lpsav = lptr(lp)
1636  lptr(lp) = lph
1637  list(lph) = in2
1638  lptr(lph) = lpsav
1639 !
1640 ! Delete IO1 as a neighbor of IO2.
1641 !
1642  call lstptr ( lend(io2), in1, list, lptr, lp )
1643  lph = lptr(lp)
1644  lptr(lp) = lptr(lph)
1645 !
1646 ! If IO1 is the last neighbor of IO2, make IN1 the last neighbor.
1647 !
1648  if ( lend(io2) == lph ) then
1649  lend(io2) = lp
1650  end if
1651 !
1652 ! Insert IN1 as a neighbor of IN2 following IO2.
1653 !
1654  call lstptr ( lend(in2), io2, list, lptr, lp )
1655  lpsav = lptr(lp)
1656  lptr(lp) = lph
1657  list(lph) = in1
1658  lptr(lph) = lpsav
1659  lp21 = lph
1660 
1661  return
1662 end subroutine swap
1663 subroutine swptst ( n1, n2, n3, n4, x, y, z, output )
1664 
1665 !*****************************************************************************80
1666 !
1667 ! Subroutine: swptst
1668 !> Decide whether to replace a diagonal arc by the other
1669 !
1670 ! Discussion:
1671 !
1672 ! This function decides whether or not to replace a
1673 ! diagonal arc in a quadrilateral with the other diagonal.
1674 ! The decision will be to swap (SWPTST = TRUE) if and only
1675 ! if N4 lies above the plane (in the half-space not contain-
1676 ! ing the origin) defined by (N1,N2,N3), or equivalently, if
1677 ! the projection of N4 onto this plane is interior to the
1678 ! circumcircle of (N1,N2,N3). The decision will be for no
1679 ! swap if the quadrilateral is not strictly convex.
1680 !
1681 ! Modified:
1682 !
1683 ! 16 June 2007
1684 !
1685 ! Author:
1686 !
1687 ! Robert Renka
1688 !
1689 ! Reference:
1690 !
1691 ! Robert Renka,
1692 ! Algorithm 772: STRIPACK,
1693 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
1694 ! ACM Transactions on Mathematical Software,
1695 ! Volume 23, Number 3, September 1997, pages 416-434.
1696 !
1697 ! Parameters:
1698 !
1699 ! Input, integer N1, N2, N3, N4, indexes of the four nodes
1700 ! defining the quadrilateral with N1 adjacent to N2, and (N1,N2,N3) in
1701 ! counterclockwise order. The arc connecting N1 to N2 should be replaced
1702 ! by an arc connecting N3 to N4 if SWPTST = TRUE. Refer to subroutine SWAP.
1703 !
1704 ! Input, real ( kind_real ) X(N), Y(N), Z(N), the coordinates of the nodes.
1705 !
1706 ! Output, logical SWPTST, TRUE if and only if the arc connecting N1
1707 ! and N2 should be swapped for an arc connecting N3 and N4.
1708 !
1709 ! Local parameters:
1710 !
1711 ! DX1,DY1,DZ1 = Coordinates of N4->N1
1712 ! DX2,DY2,DZ2 = Coordinates of N4->N2
1713 ! DX3,DY3,DZ3 = Coordinates of N4->N3
1714 ! X4,Y4,Z4 = Coordinates of N4
1715 !
1716  implicit none
1717 
1718  real ( kind_real ) dx1
1719  real ( kind_real ) dx2
1720  real ( kind_real ) dx3
1721  real ( kind_real ) dy1
1722  real ( kind_real ) dy2
1723  real ( kind_real ) dy3
1724  real ( kind_real ) dz1
1725  real ( kind_real ) dz2
1726  real ( kind_real ) dz3
1727  integer n1
1728  integer n2
1729  integer n3
1730  integer n4
1731  logical output
1732  real ( kind_real ) x(*)
1733  real ( kind_real ) x4
1734  real ( kind_real ) y(*)
1735  real ( kind_real ) y4
1736  real ( kind_real ) z(*)
1737  real ( kind_real ) z4
1738  real ( kind_real ) zz
1739 
1740  x4 = x(n4)
1741  y4 = y(n4)
1742  z4 = z(n4)
1743  dx1 = x(n1) - x4
1744  dx2 = x(n2) - x4
1745  dx3 = x(n3) - x4
1746  dy1 = y(n1) - y4
1747  dy2 = y(n2) - y4
1748  dy3 = y(n3) - y4
1749  dz1 = z(n1) - z4
1750  dz2 = z(n2) - z4
1751  dz3 = z(n3) - z4
1752 !
1753 ! N4 lies above the plane of (N1,N2,N3) iff N3 lies above
1754 ! the plane of (N2,N1,N4) iff Det(N3-N4,N2-N4,N1-N4) =
1755 ! (N3-N4,N2-N4 X N1-N4) > 0.
1756 !
1757 
1758  call det ( dx2, dy2, dz2, dx1, dy1, dz1, dx3, dy3, dz3, zz )
1759  output = sup(zz,0.0_kind_real)
1760 
1761  return
1762 end subroutine swptst
1763 subroutine trfind ( nst, p, n, x, y, z, list, lptr, lend, b1, b2, b3, i1, &
1764  i2, i3 )
1765 
1766 !*****************************************************************************80
1767 !
1768 ! Subroutine: trfind
1769 !> Locate a point relative to a triangulation
1770 !
1771 ! Discussion:
1772 !
1773 ! This subroutine locates a point P relative to a triangulation
1774 ! created by TRMESH. If P is contained in
1775 ! a triangle, the three vertex indexes and barycentric
1776 ! coordinates are returned. Otherwise, the indexes of the
1777 ! visible boundary nodes are returned.
1778 !
1779 ! Modified:
1780 !
1781 ! 16 June 2007
1782 !
1783 ! Author:
1784 !
1785 ! Robert Renka
1786 !
1787 ! Reference:
1788 !
1789 ! Robert Renka,
1790 ! Algorithm 772: STRIPACK,
1791 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
1792 ! ACM Transactions on Mathematical Software,
1793 ! Volume 23, Number 3, September 1997, pages 416-434.
1794 !
1795 ! Parameters:
1796 !
1797 ! Input, integer NST, index of a node at which TRFIND begins
1798 ! its search. Search time depends on the proximity of this node to P.
1799 !
1800 ! Input, real ( kind_real ) P(3), the x, y, and z coordinates (in that order)
1801 ! of the point P to be located.
1802 !
1803 ! Input, integer N, the number of nodes in the triangulation.
1804 ! 3 <= N.
1805 !
1806 ! Input, real ( kind_real ) X(N), Y(N), Z(N), the coordinates of the
1807 ! triangulation nodes (unit vectors).
1808 !
1809 ! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), the
1810 ! data structure defining the triangulation, created by TRMESH.
1811 !
1812 ! Output, real ( kind_real ) B1, B2, B3, the unnormalized barycentric
1813 ! coordinates of the central projection of P onto the underlying planar
1814 ! triangle if P is in the convex hull of the nodes. These parameters
1815 ! are not altered if I1 = 0.
1816 !
1817 ! Output, integer I1, I2, I3, the counterclockwise-ordered
1818 ! vertex indexes of a triangle containing P if P is contained in a triangle.
1819 ! If P is not in the convex hull of the nodes, I1 and I2 are the rightmost
1820 ! and leftmost (boundary) nodes that are visible from P, and I3 = 0. (If
1821 ! all boundary nodes are visible from P, then I1 and I2 coincide.)
1822 ! I1 = I2 = I3 = 0 if P and all of the nodes are coplanar (lie on a
1823 ! common great circle.
1824 !
1825 ! Local parameters:
1826 !
1827 ! EPS = Machine precision
1828 ! IX,IY,IZ = Integer seeds for JRAND
1829 ! LP = LIST pointer
1830 ! N0,N1,N2 = Nodes in counterclockwise order defining a
1831 ! cone (with vertex N0) containing P, or end-
1832 ! points of a boundary edge such that P Right
1833 ! N1->N2
1834 ! N1S,N2S = Initially-determined values of N1 and N2
1835 ! N3,N4 = Nodes opposite N1->N2 and N2->N1, respectively
1836 ! NEXT = Candidate for I1 or I2 when P is exterior
1837 ! NF,NL = First and last neighbors of N0, or first
1838 ! (rightmost) and last (leftmost) nodes
1839 ! visible from P when P is exterior to the
1840 ! triangulation
1841 ! PTN1 = Scalar product <P,N1>
1842 ! PTN2 = Scalar product <P,N2>
1843 ! Q = (N2 X N1) X N2 or N1 X (N2 X N1) -- used in
1844 ! the boundary traversal when P is exterior
1845 ! S12 = Scalar product <N1,N2>
1846 ! TOL = Tolerance (multiple of EPS) defining an upper
1847 ! bound on the magnitude of a negative bary-
1848 ! centric coordinate (B1 or B2) for P in a
1849 ! triangle -- used to avoid an infinite number
1850 ! of restarts with 0 <= B3 < EPS and B1 < 0 or
1851 ! B2 < 0 but small in magnitude
1852 ! XP,YP,ZP = Local variables containing P(1), P(2), and P(3)
1853 ! X0,Y0,Z0 = Dummy arguments for DET
1854 ! X1,Y1,Z1 = Dummy arguments for DET
1855 ! X2,Y2,Z2 = Dummy arguments for DET
1856 !
1857  implicit none
1858 
1859  integer n
1860 
1861  real ( kind_real ) b1
1862  real ( kind_real ) b2
1863  real ( kind_real ) b3
1864  real ( kind_real ) output
1865  real ( kind_real ) eps
1866  integer i1
1867  integer i2
1868  integer i3
1869  integer ( kind = 4 ), save :: ix = 1
1870  integer ( kind = 4 ), save :: iy = 2
1871  integer ( kind = 4 ), save :: iz = 3
1872  integer lend(n)
1873  integer list(6*(n-2))
1874  integer lp
1875  integer lptr(6*(n-2))
1876  integer n0
1877  integer n1
1878  integer n1s
1879  integer n2
1880  integer n2s
1881  integer n3
1882  integer n4
1883  integer next
1884  integer nf
1885  integer nl
1886  integer nst
1887  real ( kind_real ) p(3)
1888  real ( kind_real ) ptn1
1889  real ( kind_real ) ptn2
1890  real ( kind_real ) q(3)
1891  real ( kind_real ) s12
1892  real ( kind_real ) tol
1893  real ( kind_real ) x(n)
1894  real ( kind_real ) xp
1895  real ( kind_real ) y(n)
1896  real ( kind_real ) yp
1897  real ( kind_real ) z(n)
1898  real ( kind_real ) zp
1899 !
1900 ! Initialize variables.
1901 !
1902  xp = p(1)
1903  yp = p(2)
1904  zp = p(3)
1905  n0 = nst
1906 
1907  if ( n0 < 1 .or. n < n0 ) then
1908  call jrand ( n, ix, iy, iz, n0 )
1909  end if
1910 !
1911 ! Compute the relative machine precision EPS and TOL.
1912 !
1913  eps = epsilon( eps )
1914  tol = 100.0_kind_real * eps
1915 !
1916 ! Set NF and NL to the first and last neighbors of N0, and initialize N1 = NF.
1917 !
1918 2 continue
1919 
1920  lp = lend(n0)
1921  nl = list(lp)
1922  lp = lptr(lp)
1923  nf = list(lp)
1924  n1 = nf
1925 !
1926 ! Find a pair of adjacent neighbors N1,N2 of N0 that define
1927 ! a wedge containing P: P LEFT N0->N1 and P RIGHT N0->N2.
1928 !
1929  if ( 0 < nl ) then
1930 !
1931 ! N0 is an interior node. Find N1.
1932 !
1933 3 continue
1934 
1935  call det(x(n0),y(n0),z(n0),x(n1),y(n1),z(n1),xp,yp,zp,output)
1936  if ( inf(output,0.0_kind_real) ) then
1937  lp = lptr(lp)
1938  n1 = list(lp)
1939  if ( n1 == nl ) then
1940  go to 6
1941  end if
1942  go to 3
1943  end if
1944  else
1945 !
1946 ! N0 is a boundary node. Test for P exterior.
1947 !
1948  nl = -nl
1949 !
1950 ! Is P to the right of the boundary edge N0->NF?
1951 !
1952  call det(x(n0),y(n0),z(n0),x(nf),y(nf),z(nf), xp,yp,zp,output)
1953  if ( inf(output,0.0_kind_real) ) then
1954  n1 = n0
1955  n2 = nf
1956  go to 9
1957  end if
1958 !
1959 ! Is P to the right of the boundary edge NL->N0?
1960 !
1961  call det(x(nl),y(nl),z(nl),x(n0),y(n0),z(n0),xp,yp,zp,output)
1962  if ( inf(output,0.0_kind_real) ) then
1963  n1 = nl
1964  n2 = n0
1965  go to 9
1966  end if
1967 
1968  end if
1969 !
1970 ! P is to the left of arcs N0->N1 and NL->N0. Set N2 to the
1971 ! next neighbor of N0 (following N1).
1972 !
1973 4 continue
1974 
1975  lp = lptr(lp)
1976  n2 = abs( list(lp) )
1977 
1978  call det(x(n0),y(n0),z(n0),x(n2),y(n2),z(n2),xp,yp,zp,output)
1979  if ( inf(output,0.0_kind_real) ) then
1980  go to 7
1981  end if
1982 
1983  n1 = n2
1984 
1985  if ( n1 /= nl ) then
1986  go to 4
1987  end if
1988 
1989  call det ( x(n0), y(n0), z(n0), x(nf), y(nf), z(nf), xp, yp, zp, output )
1990  if ( inf(output,0.0_kind_real) ) then
1991  go to 6
1992  end if
1993 !
1994 ! P is left of or on arcs N0->NB for all neighbors NB
1995 ! of N0. Test for P = +/-N0.
1996 !
1997  if ( inf(abs( x(n0 ) * xp + y(n0) * yp + z(n0) * zp),1.0-4.0*eps) ) then
1998 !
1999 ! All points are collinear iff P Left NB->N0 for all
2000 ! neighbors NB of N0. Search the neighbors of N0.
2001 ! Note: N1 = NL and LP points to NL.
2002 !
2003  do
2004 
2005  call det(x(n1),y(n1),z(n1),x(n0),y(n0),z(n0),xp,yp,zp,output)
2006  if ( inf(output,0.0_kind_real) ) then
2007  exit
2008  end if
2009 
2010  lp = lptr(lp)
2011  n1 = abs( list(lp) )
2012 
2013  if ( n1 == nl ) then
2014  i1 = 0
2015  i2 = 0
2016  i3 = 0
2017  return
2018  end if
2019 
2020  end do
2021 
2022  end if
2023 !
2024 ! P is to the right of N1->N0, or P = +/-N0. Set N0 to N1 and start over.
2025 !
2026  n0 = n1
2027  go to 2
2028 !
2029 ! P is between arcs N0->N1 and N0->NF.
2030 !
2031 6 continue
2032 
2033  n2 = nf
2034 !
2035 ! P is contained in a wedge defined by geodesics N0-N1 and
2036 ! N0-N2, where N1 is adjacent to N2. Save N1 and N2 to
2037 ! test for cycling.
2038 !
2039 7 continue
2040 
2041  n3 = n0
2042  n1s = n1
2043  n2s = n2
2044 !
2045 ! Top of edge-hopping loop:
2046 !
2047 8 continue
2048 
2049  call det ( x(n1),y(n1),z(n1),x(n2),y(n2),z(n2),xp,yp,zp,b3 )
2050 
2051  if ( inf(b3,0.0_kind_real) ) then
2052 !
2053 ! Set N4 to the first neighbor of N2 following N1 (the
2054 ! node opposite N2->N1) unless N1->N2 is a boundary arc.
2055 !
2056  call lstptr ( lend(n2), n1, list, lptr, lp )
2057 
2058  if ( list(lp) < 0 ) then
2059  go to 9
2060  end if
2061 
2062  lp = lptr(lp)
2063  n4 = abs( list(lp) )
2064 !
2065 ! Define a new arc N1->N2 which intersects the geodesic N0-P.
2066 !
2067  call det ( x(n0),y(n0),z(n0),x(n4),y(n4),z(n4),xp,yp,zp,output )
2068  if ( inf(output,0.0_kind_real) ) then
2069  n3 = n2
2070  n2 = n4
2071  n1s = n1
2072  if ( n2 /= n2s .and. n2 /= n0 ) then
2073  go to 8
2074  end if
2075  else
2076  n3 = n1
2077  n1 = n4
2078  n2s = n2
2079  if ( n1 /= n1s .and. n1 /= n0 ) then
2080  go to 8
2081  end if
2082  end if
2083 !
2084 ! The starting node N0 or edge N1-N2 was encountered
2085 ! again, implying a cycle (infinite loop). Restart
2086 ! with N0 randomly selected.
2087 !
2088  call jrand ( n, ix, iy, iz, n0 )
2089  go to 2
2090 
2091  end if
2092 !
2093 ! P is in (N1,N2,N3) unless N0, N1, N2, and P are collinear
2094 ! or P is close to -N0.
2095 !
2096  if ( supeq(b3,eps) ) then
2097 !
2098 ! B3 /= 0.
2099 !
2100  call det(x(n2),y(n2),z(n2),x(n3),y(n3),z(n3),xp,yp,zp,b1)
2101  call det(x(n3),y(n3),z(n3),x(n1),y(n1),z(n1),xp,yp,zp,b2)
2102 !
2103 ! Restart with N0 randomly selected.
2104 !
2105  if ( inf(b1,-tol) .or. inf(b2,-tol) ) then
2106  call jrand ( n, ix, iy, iz, n0 )
2107  go to 2
2108  end if
2109 
2110  else
2111 !
2112 ! B3 = 0 and thus P lies on N1->N2. Compute
2113 ! B1 = Det(P,N2 X N1,N2) and B2 = Det(P,N1,N2 X N1).
2114 !
2115  b3 = 0.0_kind_real
2116  s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
2117  ptn1 = xp * x(n1) + yp * y(n1) + zp * z(n1)
2118  ptn2 = xp * x(n2) + yp * y(n2) + zp * z(n2)
2119  b1 = ptn1 - s12 * ptn2
2120  b2 = ptn2 - s12 * ptn1
2121 !
2122 ! Restart with N0 randomly selected.
2123 !
2124  if ( inf(b1,-tol) .or. inf(b2,-tol) ) then
2125  call jrand ( n, ix, iy, iz, n0 )
2126  go to 2
2127  end if
2128 
2129  end if
2130 !
2131 ! P is in (N1,N2,N3).
2132 !
2133  i1 = n1
2134  i2 = n2
2135  i3 = n3
2136  b1 = max( b1, 0.0_kind_real )
2137  b2 = max( b2, 0.0_kind_real )
2138  return
2139 !
2140 ! P Right N1->N2, where N1->N2 is a boundary edge.
2141 ! Save N1 and N2, and set NL = 0 to indicate that
2142 ! NL has not yet been found.
2143 !
2144 9 continue
2145 
2146  n1s = n1
2147  n2s = n2
2148  nl = 0
2149 !
2150 ! Counterclockwise Boundary Traversal:
2151 !
2152 10 continue
2153 
2154  lp = lend(n2)
2155  lp = lptr(lp)
2156  next = list(lp)
2157 
2158  call det(x(n2),y(n2),z(n2),x(next),y(next),z(next),xp,yp,zp,output)
2159  if ( supeq(output,0.0_kind_real) ) then
2160 !
2161 ! N2 is the rightmost visible node if P Forward N2->N1
2162 ! or NEXT Forward N2->N1. Set Q to (N2 X N1) X N2.
2163 !
2164  s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
2165 
2166  q(1) = x(n1) - s12 * x(n2)
2167  q(2) = y(n1) - s12 * y(n2)
2168  q(3) = z(n1) - s12 * z(n2)
2169 
2170  if ( sup(xp * q(1) + yp * q(2) + zp * q(3),0.0_kind_real) ) then
2171  go to 11
2172  end if
2173 
2174  if ( sup(x(next) * q(1) + y(next) * q(2) + z(next) * q(3),0.0_kind_real) ) then
2175  go to 11
2176  end if
2177 !
2178 ! N1, N2, NEXT, and P are nearly collinear, and N2 is
2179 ! the leftmost visible node.
2180 !
2181  nl = n2
2182  end if
2183 !
2184 ! Bottom of counterclockwise loop:
2185 !
2186  n1 = n2
2187  n2 = next
2188 
2189  if ( n2 /= n1s ) then
2190  go to 10
2191  end if
2192 !
2193 ! All boundary nodes are visible from P.
2194 !
2195  i1 = n1s
2196  i2 = n1s
2197  i3 = 0
2198  return
2199 !
2200 ! N2 is the rightmost visible node.
2201 !
2202 11 continue
2203 
2204  nf = n2
2205 
2206  if ( nl == 0 ) then
2207 !
2208 ! Restore initial values of N1 and N2, and begin the search
2209 ! for the leftmost visible node.
2210 !
2211  n2 = n2s
2212  n1 = n1s
2213 !
2214 ! Clockwise Boundary Traversal:
2215 !
2216 12 continue
2217 
2218  lp = lend(n1)
2219  next = -list(lp)
2220 
2221  call det ( x(next), y(next), z(next), x(n1), y(n1), z(n1), xp, yp, zp, output )
2222  if ( supeq(output,0.0_kind_real) ) then
2223 !
2224 ! N1 is the leftmost visible node if P or NEXT is
2225 ! forward of N1->N2. Compute Q = N1 X (N2 X N1).
2226 !
2227  s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
2228  q(1) = x(n2) - s12 * x(n1)
2229  q(2) = y(n2) - s12 * y(n1)
2230  q(3) = z(n2) - s12 * z(n1)
2231 
2232  if ( sup(xp * q(1) + yp * q(2) + zp * q(3),0.0_kind_real) ) then
2233  go to 13
2234  end if
2235 
2236  if ( sup(x(next) * q(1) + y(next) * q(2) + z(next) * q(3),0.0_kind_real) ) then
2237  go to 13
2238  end if
2239 !
2240 ! P, NEXT, N1, and N2 are nearly collinear and N1 is the rightmost
2241 ! visible node.
2242 !
2243  nf = n1
2244  end if
2245 !
2246 ! Bottom of clockwise loop:
2247 !
2248  n2 = n1
2249  n1 = next
2250 
2251  if ( n1 /= n1s ) then
2252  go to 12
2253  end if
2254 !
2255 ! All boundary nodes are visible from P.
2256 !
2257  i1 = n1
2258  i2 = n1
2259  i3 = 0
2260  return
2261 !
2262 ! N1 is the leftmost visible node.
2263 !
2264 13 continue
2265 
2266  nl = n1
2267 
2268  end if
2269 !
2270 ! NF and NL have been found.
2271 !
2272  i1 = nf
2273  i2 = nl
2274  i3 = 0
2275 
2276  return
2277 end subroutine trfind
2278 subroutine trlist ( n, list, lptr, lend, nrow, nt, ltri, ier )
2279 
2280 !*****************************************************************************80
2281 !
2282 ! Subroutine: trlist
2283 !> Convert a triangulation data structure to a triangle list
2284 !
2285 ! Discussion:
2286 !
2287 ! This subroutine converts a triangulation data structure
2288 ! from the linked list created by TRMESH to a triangle list.
2289 !
2290 ! Modified:
2291 !
2292 ! 16 June 2007
2293 !
2294 ! Author:
2295 !
2296 ! Robert Renka
2297 !
2298 ! Reference:
2299 !
2300 ! Robert Renka,
2301 ! Algorithm 772: STRIPACK,
2302 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
2303 ! ACM Transactions on Mathematical Software,
2304 ! Volume 23, Number 3, September 1997, pages 416-434.
2305 !
2306 ! Parameters:
2307 !
2308 ! Input, integer N, the number of nodes in the triangulation.
2309 ! 3 <= N.
2310 !
2311 ! Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), linked
2312 ! list data structure defining the triangulation. Refer to TRMESH.
2313 !
2314 ! Input, integer NROW, the number of rows (entries per triangle)
2315 ! reserved for the triangle list LTRI. The value must be 6 if only the
2316 ! vertex indexes and neighboring triangle indexes are to be stored, or 9
2317 ! if arc indexes are also to be assigned and stored. Refer to LTRI.
2318 !
2319 ! Output, integer NT, the number of triangles in the
2320 ! triangulation unless IER /=0, in which case NT = 0. NT = 2N-NB-2 if
2321 ! NB >= 3 or 2N-4 if NB = 0, where NB is the number of boundary nodes.
2322 !
2323 ! Output, integer LTRI(NROW,*). The second dimension of LTRI
2324 ! must be at least NT, where NT will be at most 2*N-4. The J-th column
2325 ! contains the vertex nodal indexes (first three rows), neighboring triangle
2326 ! indexes (second three rows), and, if NROW = 9, arc indexes (last three
2327 ! rows) associated with triangle J for J = 1,...,NT. The vertices are
2328 ! ordered counterclockwise with the first vertex taken to be the one
2329 ! with smallest index. Thus, LTRI(2,J) and LTRI(3,J) are larger than
2330 ! LTRI(1,J) and index adjacent neighbors of node LTRI(1,J). For
2331 ! I = 1,2,3, LTRI(I+3,J) and LTRI(I+6,J) index the triangle and arc,
2332 ! respectively, which are opposite (not shared by) node LTRI(I,J), with
2333 ! LTRI(I+3,J) = 0 if LTRI(I+6,J) indexes a boundary arc. Vertex indexes
2334 ! range from 1 to N, triangle indexes from 0 to NT, and, if included,
2335 ! arc indexes from 1 to NA, where NA = 3N-NB-3 if NB >= 3 or 3N-6 if
2336 ! NB = 0. The triangles are ordered on first (smallest) vertex indexes.
2337 !
2338 ! Output, integer IER, error indicator.
2339 ! 0, if no errors were encountered.
2340 ! 1, if N or NROW is outside its valid range on input.
2341 ! 2, if the triangulation data structure (LIST,LPTR,LEND) is invalid.
2342 ! Note, however, that these arrays are not completely tested for validity.
2343 !
2344 ! Local parameters:
2345 !
2346 ! ARCS = Logical variable with value TRUE iff are
2347 ! indexes are to be stored
2348 ! I,J = LTRI row indexes (1 to 3) associated with
2349 ! triangles KT and KN, respectively
2350 ! I1,I2,I3 = Nodal indexes of triangle KN
2351 ! ISV = Variable used to permute indexes I1,I2,I3
2352 ! KA = Arc index and number of currently stored arcs
2353 ! KN = Index of the triangle that shares arc I1-I2 with KT
2354 ! KT = Triangle index and number of currently stored triangles
2355 ! LP = LIST pointer
2356 ! LP2 = Pointer to N2 as a neighbor of N1
2357 ! LPL = Pointer to the last neighbor of I1
2358 ! LPLN1 = Pointer to the last neighbor of N1
2359 ! N1,N2,N3 = Nodal indexes of triangle KT
2360 ! NM2 = N-2
2361 !
2362  implicit none
2363 
2364  integer n
2365  integer nrow
2366 
2367  logical arcs
2368  integer i
2369  integer i1
2370  integer i2
2371  integer i3
2372  integer ier
2373  integer isv
2374  integer j
2375  integer ka
2376  integer kn
2377  integer kt
2378  integer lend(n)
2379  integer list(6*(n-2))
2380  integer lp
2381  integer lp2
2382  integer lpl
2383  integer lpln1
2384  integer lptr(6*(n-2))
2385  integer ltri(nrow,*)
2386  integer n1
2387  integer n2
2388  integer n3
2389  integer nm2
2390  integer nt
2391 !
2392 ! Test for invalid input parameters.
2393 !
2394  if ( n < 3 .or. ( nrow /= 6 .and. nrow /= 9 ) ) then
2395  nt = 0
2396  ier = 1
2397  return
2398  end if
2399 !
2400 ! Initialize parameters for loop on triangles KT = (N1,N2,
2401 ! N3), where N1 < N2 and N1 < N3.
2402 !
2403 ! ARCS = TRUE iff arc indexes are to be stored.
2404 ! KA,KT = Numbers of currently stored arcs and triangles.
2405 ! NM2 = Upper bound on candidates for N1.
2406 !
2407  arcs = nrow == 9
2408  ka = 0
2409  kt = 0
2410  nm2 = n-2
2411 !
2412 ! Loop on nodes N1.
2413 !
2414  do n1 = 1, nm2
2415 !
2416 ! Loop on pairs of adjacent neighbors (N2,N3). LPLN1 points
2417 ! to the last neighbor of N1, and LP2 points to N2.
2418 !
2419  lpln1 = lend(n1)
2420  lp2 = lpln1
2421 
2422 1 continue
2423 
2424  lp2 = lptr(lp2)
2425  n2 = list(lp2)
2426  lp = lptr(lp2)
2427  n3 = abs( list(lp) )
2428 
2429  if ( n2 < n1 .or. n3 < n1 ) then
2430  go to 8
2431  end if
2432 !
2433 ! Add a new triangle KT = (N1,N2,N3).
2434 !
2435  kt = kt + 1
2436  ltri(1,kt) = n1
2437  ltri(2,kt) = n2
2438  ltri(3,kt) = n3
2439 !
2440 ! Loop on triangle sides (I2,I1) with neighboring triangles
2441 ! KN = (I1,I2,I3).
2442 !
2443  do i = 1, 3
2444 
2445  if ( i == 1 ) then
2446  i1 = n3
2447  i2 = n2
2448  else if ( i == 2 ) then
2449  i1 = n1
2450  i2 = n3
2451  else
2452  i1 = n2
2453  i2 = n1
2454  end if
2455  j = 0
2456 !
2457 ! Set I3 to the neighbor of I1 that follows I2 unless
2458 ! I2->I1 is a boundary arc.
2459 !
2460  lpl = lend(i1)
2461  lp = lptr(lpl)
2462 
2463  do
2464 
2465  if ( list(lp) == i2 ) then
2466  go to 3
2467  end if
2468 
2469  lp = lptr(lp)
2470 
2471  if ( lp == lpl ) then
2472  exit
2473  end if
2474 
2475  end do
2476 !
2477 ! Invalid triangulation data structure: I1 is a neighbor of
2478 ! I2, but I2 is not a neighbor of I1.
2479 !
2480  if ( abs( list(lp) ) /= i2 ) then
2481  nt = 0
2482  ier = 2
2483  return
2484  end if
2485 !
2486 ! I2 is the last neighbor of I1. Bypass the search for a neighboring
2487 ! triangle if I2->I1 is a boundary arc.
2488 !
2489  kn = 0
2490 
2491  if ( list(lp) < 0 ) then
2492  go to 6
2493  end if
2494 !
2495 ! I2->I1 is not a boundary arc, and LP points to I2 as
2496 ! a neighbor of I1.
2497 !
2498 3 continue
2499 
2500  lp = lptr(lp)
2501  i3 = abs( list(lp) )
2502 !
2503 ! Find J such that LTRI(J,KN) = I3 (not used if KT < KN),
2504 ! and permute the vertex indexes of KN so that I1 is smallest.
2505 !
2506  if ( i1 < i2 .and. i1 < i3 ) then
2507  j = 3
2508  else if ( i2 < i3 ) then
2509  j = 2
2510  isv = i1
2511  i1 = i2
2512  i2 = i3
2513  i3 = isv
2514  else
2515  j = 1
2516  isv = i1
2517  i1 = i3
2518  i3 = i2
2519  i2 = isv
2520  end if
2521 !
2522 ! Test for KT < KN (triangle index not yet assigned).
2523 !
2524  if ( n1 < i1 ) then
2525  cycle
2526  end if
2527 !
2528 ! Find KN, if it exists, by searching the triangle list in
2529 ! reverse order.
2530 !
2531  do kn = kt-1, 1, -1
2532  if ( ltri(1,kn) == i1 .and. &
2533  ltri(2,kn) == i2 .and. &
2534  ltri(3,kn) == i3 ) then
2535  go to 5
2536  end if
2537  end do
2538 
2539  cycle
2540 !
2541 ! Store KT as a neighbor of KN.
2542 !
2543 5 continue
2544 
2545  ltri(j+3,kn) = kt
2546 !
2547 ! Store KN as a neighbor of KT, and add a new arc KA.
2548 !
2549 6 continue
2550 
2551  ltri(i+3,kt) = kn
2552 
2553  if ( arcs ) then
2554  ka = ka + 1
2555  ltri(i+6,kt) = ka
2556  if ( kn /= 0 ) then
2557  ltri(j+6,kn) = ka
2558  end if
2559  end if
2560 
2561  end do
2562 !
2563 ! Bottom of loop on triangles.
2564 !
2565 8 continue
2566 
2567  if ( lp2 /= lpln1 ) then
2568  go to 1
2569  end if
2570 
2571  end do
2572 
2573  nt = kt
2574  ier = 0
2575 
2576  return
2577 end subroutine trlist
2578 subroutine trmesh ( mpl, n, x, y, z, list, lptr, lend, lnew, near, next, dist, ier )
2579 
2580 !*****************************************************************************80
2581 !
2582 ! Subroutine: trmesh
2583 !> Create a Delaunay triangulation on the unit sphere
2584 !
2585 ! Discussion:
2586 !
2587 ! This subroutine creates a Delaunay triangulation of a
2588 ! set of N arbitrarily distributed points, referred to as
2589 ! nodes, on the surface of the unit sphere. The Delaunay
2590 ! triangulation is defined as a set of (spherical) triangles
2591 ! with the following five properties:
2592 !
2593 ! 1) The triangle vertices are nodes.
2594 ! 2) No triangle contains a node other than its vertices.
2595 ! 3) The interiors of the triangles are pairwise disjoint.
2596 ! 4) The union of triangles is the convex hull of the set
2597 ! of nodes (the smallest convex set that contains
2598 ! the nodes). If the nodes are not contained in a
2599 ! single hemisphere, their convex hull is the
2600 ! entire sphere and there are no boundary nodes.
2601 ! Otherwise, there are at least three boundary nodes.
2602 ! 5) The interior of the circumcircle of each triangle
2603 ! contains no node.
2604 !
2605 ! The first four properties define a triangulation, and the
2606 ! last property results in a triangulation which is as close
2607 ! as possible to equiangular in a certain sense and which is
2608 ! uniquely defined unless four or more nodes lie in a common
2609 ! plane. This property makes the triangulation well-suited
2610 ! for solving closest-point problems and for triangle-based
2611 ! interpolation.
2612 !
2613 ! Provided the nodes are randomly ordered, the algorithm
2614 ! has expected time complexity O(N*log(N)) for most nodal
2615 ! distributions. Note, however, that the complexity may be
2616 ! as high as O(N**2) if, for example, the nodes are ordered
2617 ! on increasing latitude.
2618 !
2619 ! Spherical coordinates (latitude and longitude) may be
2620 ! converted to Cartesian coordinates by Subroutine TRANS.
2621 !
2622 ! The following is a list of the software package modules
2623 ! which a user may wish to call directly:
2624 !
2625 ! ADDNOD - Updates the triangulation by appending a new node.
2626 !
2627 ! BNODES - Returns an array containing the indexes of the
2628 ! boundary nodes (if any) in counterclockwise
2629 ! order. Counts of boundary nodes, triangles,
2630 ! and arcs are also returned.
2631 !
2632 ! EDGE - Forces an arbitrary pair of nodes to be connected
2633 ! by an arc in the triangulation.
2634 !
2635 ! INSIDE - Locates a point relative to a polygon on the
2636 ! surface of the sphere.
2637 !
2638 ! INTRSC - Returns the point of intersection between a
2639 ! pair of great circle arcs.
2640 !
2641 ! JRAND - Generates a uniformly distributed pseudo-random integer.
2642 !
2643 ! LEFT - Locates a point relative to a great circle.
2644 !
2645 ! TRANS - Transforms spherical coordinates into Cartesian
2646 ! coordinates on the unit sphere for input to
2647 ! Subroutine TRMESH.
2648 !
2649 ! TRLIST - Converts the triangulation data structure to a
2650 ! triangle list more suitable for use in a finite
2651 ! element code.
2652 !
2653 ! TRMESH - Creates a Delaunay triangulation of a set of
2654 ! nodes.
2655 !
2656 ! Modified:
2657 !
2658 ! 16 June 2007
2659 !
2660 ! Author:
2661 !
2662 ! Robert Renka
2663 !
2664 ! Reference:
2665 !
2666 ! Robert Renka,
2667 ! Algorithm 772: STRIPACK,
2668 ! Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
2669 ! ACM Transactions on Mathematical Software,
2670 ! Volume 23, Number 3, September 1997, pages 416-434.
2671 !
2672 ! Parameters:
2673 !
2674 ! Input, integer N, the number of nodes in the triangulation.
2675 ! 3 <= N.
2676 !
2677 ! Input, real ( kind_real ) X(N), Y(N), Z(N), the coordinates of distinct
2678 ! nodes. (X(K),Y(K), Z(K)) is referred to as node K, and K is referred
2679 ! to as a nodal index. It is required that X(K)**2 + Y(K)**2 + Z(K)**2 = 1
2680 ! for all K. The first three nodes must not be collinear (lie on a
2681 ! common great circle).
2682 !
2683 ! Output, integer LIST(6*(N-2)), nodal indexes which, along
2684 ! with LPTR, LEND, and LNEW, define the triangulation as a set of N
2685 ! adjacency lists; counterclockwise-ordered sequences of neighboring nodes
2686 ! such that the first and last neighbors of a boundary node are boundary
2687 ! nodes (the first neighbor of an interior node is arbitrary). In order to
2688 ! distinguish between interior and boundary nodes, the last neighbor of
2689 ! each boundary node is represented by the negative of its index.
2690 !
2691 ! Output, integer LPTR(6*(N-2)), = Set of pointers (LIST
2692 ! indexes) in one-to-one correspondence with the elements of LIST.
2693 ! LIST(LPTR(I)) indexes the node which follows LIST(I) in cyclical
2694 ! counterclockwise order (the first neighbor follows the last neighbor).
2695 !
2696 ! Output, integer LEND(N), pointers to adjacency lists.
2697 ! LEND(K) points to the last neighbor of node K. LIST(LEND(K)) < 0 if and
2698 ! only if K is a boundary node.
2699 !
2700 ! Output, integer LNEW, pointer to the first empty location
2701 ! in LIST and LPTR (list length plus one). LIST, LPTR, LEND, and LNEW are
2702 ! not altered if IER < 0, and are incomplete if 0 < IER.
2703 !
2704 ! Workspace, integer NEAR(N),
2705 ! used to efficiently determine the nearest triangulation node to each
2706 ! unprocessed node for use by ADDNOD.
2707 !
2708 ! Workspace, integer NEXT(N),
2709 ! used to efficiently determine the nearest triangulation node to each
2710 ! unprocessed node for use by ADDNOD.
2711 !
2712 ! Workspace, real ( kind_real ) DIST(N),
2713 ! used to efficiently determine the nearest triangulation node to each
2714 ! unprocessed node for use by ADDNOD.
2715 !
2716 ! Output, integer IER, error indicator:
2717 ! 0, if no errors were encountered.
2718 ! -1, if N < 3 on input.
2719 ! -2, if the first three nodes are collinear.
2720 ! L, if nodes L and M coincide for some L < M. The data structure
2721 ! represents a triangulation of nodes 1 to M-1 in this case.
2722 !
2723 ! Local parameters:
2724 !
2725 ! D = (Negative cosine of) distance from node K to node I
2726 ! D1,D2,D3 = Distances from node K to nodes 1, 2, and 3, respectively
2727 ! I,J = Nodal indexes
2728 ! I0 = Index of the node preceding I in a sequence of
2729 ! unprocessed nodes: I = NEXT(I0)
2730 ! K = Index of node to be added and DO-loop index: 3 < K
2731 ! LP = LIST index (pointer) of a neighbor of K
2732 ! LPL = Pointer to the last neighbor of K
2733 ! NEXTI = NEXT(I)
2734 ! NN = Local copy of N
2735 !
2736  implicit none
2737 
2738  type(mpl_type),intent(inout) :: mpl !< MPI data
2739 
2740  integer n
2741 
2742  real ( kind_real ) d
2743  real ( kind_real ) d1
2744  real ( kind_real ) d2
2745  real ( kind_real ) d3
2746  real ( kind_real ) dist(n)
2747  integer i
2748  integer i0
2749  integer ier
2750  integer j
2751  integer k
2752  logical output_1, output_2
2753  integer lend(n)
2754  integer list(6*(n-2))
2755  integer lnew
2756  integer lp
2757  integer lpl
2758  integer lptr(6*(n-2))
2759  integer near(n)
2760  integer next(n)
2761  integer nexti
2762  integer nn
2763  real ( kind_real ) x(n)
2764  real ( kind_real ) y(n)
2765  real ( kind_real ) z(n)
2766  character(len=1024),parameter :: subr = 'trmesh'
2767 
2768  nn = n
2769 
2770  if ( nn < 3 ) then
2771  ier = -1
2772  write ( *, '(a)' ) ' '
2773  write ( *, '(a)' ) 'TRMESH - Fatal error!'
2774  write ( *, '(a)' ) ' N < 3.'
2775  call mpl%abort(subr,'stop in stripack')
2776  end if
2777 !
2778 ! Initialize
2779 !
2780  list = 0
2781  lptr = 0
2782  lend = 0
2783 !
2784 ! Store the first triangle in the linked list.
2785 !
2786  call left (x(1),y(1),z(1),x(2),y(2),z(2),x(3),y(3),z(3),output_1)
2787  call left (x(2),y(2),z(2),x(1),y(1),z(1),x(3),y(3),z(3),output_2)
2788 
2789  if ( .not. output_1 ) then
2790 !
2791 ! The first triangle is (3,2,1) = (2,1,3) = (1,3,2).
2792 !
2793  list(1) = 3
2794  lptr(1) = 2
2795  list(2) = -2
2796  lptr(2) = 1
2797  lend(1) = 2
2798 
2799  list(3) = 1
2800  lptr(3) = 4
2801  list(4) = -3
2802  lptr(4) = 3
2803  lend(2) = 4
2804 
2805  list(5) = 2
2806  lptr(5) = 6
2807  list(6) = -1
2808  lptr(6) = 5
2809  lend(3) = 6
2810 
2811  else if ( .not. output_2 ) then
2812 !
2813 ! The first triangle is (1,2,3): 3 Strictly Left 1->2,
2814 ! i.e., node 3 lies in the left hemisphere defined by arc 1->2.
2815 !
2816  list(1) = 2
2817  lptr(1) = 2
2818  list(2) = -3
2819  lptr(2) = 1
2820  lend(1) = 2
2821 
2822  list(3) = 3
2823  lptr(3) = 4
2824  list(4) = -1
2825  lptr(4) = 3
2826  lend(2) = 4
2827 
2828  list(5) = 1
2829  lptr(5) = 6
2830  list(6) = -2
2831  lptr(6) = 5
2832  lend(3) = 6
2833 !
2834 ! The first three nodes are collinear.
2835 !
2836  else
2837 
2838  ier = -2
2839  return
2840 
2841  end if
2842 !
2843 ! Initialize LNEW and test for N = 3.
2844 !
2845  lnew = 7
2846 
2847  if ( nn == 3 ) then
2848  ier = 0
2849  return
2850  end if
2851 !
2852 ! A nearest-node data structure (NEAR, NEXT, and DIST) is
2853 ! used to obtain an expected-time (N*log(N)) incremental
2854 ! algorithm by enabling constant search time for locating
2855 ! each new node in the triangulation.
2856 !
2857 ! For each unprocessed node K, NEAR(K) is the index of the
2858 ! triangulation node closest to K (used as the starting
2859 ! point for the search in Subroutine TRFIND) and DIST(K)
2860 ! is an increasing function of the arc length (angular
2861 ! distance) between nodes K and NEAR(K): -Cos(a) for arc
2862 ! length a.
2863 !
2864 ! Since it is necessary to efficiently find the subset of
2865 ! unprocessed nodes associated with each triangulation
2866 ! node J (those that have J as their NEAR entries), the
2867 ! subsets are stored in NEAR and NEXT as follows: for
2868 ! each node J in the triangulation, I = NEAR(J) is the
2869 ! first unprocessed node in J's set (with I = 0 if the
2870 ! set is empty), L = NEXT(I) (if 0 < I) is the second,
2871 ! NEXT(L) (if 0 < L) is the third, etc. The nodes in each
2872 ! set are initially ordered by increasing indexes (which
2873 ! maximizes efficiency) but that ordering is not main-
2874 ! tained as the data structure is updated.
2875 !
2876 ! Initialize the data structure for the single triangle.
2877 !
2878  near(1) = 0
2879  near(2) = 0
2880  near(3) = 0
2881 
2882  do k = nn, 4, -1
2883 
2884  d1 = -( x(k) * x(1) + y(k) * y(1) + z(k) * z(1) )
2885  d2 = -( x(k) * x(2) + y(k) * y(2) + z(k) * z(2) )
2886  d3 = -( x(k) * x(3) + y(k) * y(3) + z(k) * z(3) )
2887 
2888  if ( inf(d1,d2) .and. inf(d1,d3) ) then
2889  near(k) = 1
2890  dist(k) = d1
2891  next(k) = near(1)
2892  near(1) = k
2893  else if ( inf(d2,d1) .and. inf(d2,d3) ) then
2894  near(k) = 2
2895  dist(k) = d2
2896  next(k) = near(2)
2897  near(2) = k
2898  else
2899  near(k) = 3
2900  dist(k) = d3
2901  next(k) = near(3)
2902  near(3) = k
2903  end if
2904  end do
2905 !
2906 ! Add the remaining nodes.
2907 !
2908  do k = 4, nn
2909 
2910  call addnod ( mpl, near(k), k, x, y, z, list, lptr, lend, lnew, ier )
2911 
2912  if ( ier /= 0 ) then
2913  write ( *, '(a)' ) 'TRMESH - ADDNOD - Fatal error!'
2914  write ( *, '(a,i8,a,i8)' ) ' Node ', k, ' is equal to already used node ', ier
2915  call mpl%abort(subr,'stop in stripack')
2916  end if
2917 !
2918 ! Remove K from the set of unprocessed nodes associated with NEAR(K).
2919 !
2920  i = near(k)
2921  i0 = i
2922 
2923  if ( near(i) == k ) then
2924 
2925  near(i) = next(k)
2926 
2927  else
2928 
2929  i = near(i)
2930 
2931  do
2932 
2933  i0 = i
2934  i = next(i0)
2935 
2936  if ( i == k ) then
2937  exit
2938  end if
2939 
2940  end do
2941 
2942  next(i0) = next(k)
2943 
2944  end if
2945 
2946  near(k) = 0
2947 !
2948 ! Loop on neighbors J of node K.
2949 !
2950  lpl = lend(k)
2951  lp = lpl
2952 
2953 3 continue
2954 
2955  lp = lptr(lp)
2956  j = abs( list(lp) )
2957 !
2958 ! Loop on elements I in the sequence of unprocessed nodes
2959 ! associated with J: K is a candidate for replacing J
2960 ! as the nearest triangulation node to I. The next value
2961 ! of I in the sequence, NEXT(I), must be saved before I
2962 ! is moved because it is altered by adding I to K's set.
2963 !
2964  i = near(j)
2965 
2966  do
2967 
2968  if ( i == 0 ) then
2969  exit
2970  end if
2971 
2972  nexti = next(i)
2973 !
2974 ! Test for the distance from I to K less than the distance
2975 ! from I to J. Indistinguishability threshold for cross-plateform reproducibility
2976 !
2977  d = - ( x(i) * x(k) + y(i) * y(k) + z(i) * z(k) )
2978  if ( inf(d,dist(i)) ) then
2979 !
2980 ! Replace J by K as the nearest triangulation node to I:
2981 ! update NEAR(I) and DIST(I), and remove I from J's set
2982 ! of unprocessed nodes and add it to K's set.
2983 !
2984  near(i) = k
2985  dist(i) = d
2986 
2987  if ( i == near(j) ) then
2988  near(j) = nexti
2989  else
2990  next(i0) = nexti
2991  end if
2992 
2993  next(i) = near(k)
2994  near(k) = i
2995  else
2996  i0 = i
2997  end if
2998 
2999  i = nexti
3000 
3001  end do
3002 !
3003 ! Bottom of loop on neighbors J.
3004 !
3005  if ( lp /= lpl ) then
3006  go to 3
3007  end if
3008  end do
3009 
3010  return
3011 end subroutine trmesh
3012 
3013 end module tools_stripack
tools_repro::sup
logical function, public sup(x, y)
Superior test for reals.
Definition: tools_repro.F90:92
tools_stripack::swptst
subroutine swptst(n1, n2, n3, n4, x, y, z, output)
Definition: tools_stripack.F90:1664
tools_repro::inf
logical function, public inf(x, y)
Inferior test for reals.
Definition: tools_repro.F90:53
tools_stripack::inside
logical function, public inside(p, lv, xv, yv, zv, nv, listv, ier)
Definition: tools_stripack.F90:764
tools_repro::supeq
logical function, public supeq(x, y)
Superior or equal test for reals.
Definition: tools_repro.F90:112
tools_stripack::bdyadd
subroutine bdyadd(kk, i1, i2, list, lptr, lend, lnew)
Definition: tools_stripack.F90:279
tools_stripack::covsph
subroutine covsph(kk, n0, list, lptr, lend, lnew)
Definition: tools_stripack.F90:557
tools_stripack::intrsc
subroutine intrsc(p1, p2, cn, p, ier)
Definition: tools_stripack.F90:1229
tools_stripack::swap
subroutine swap(in1, in2, io1, io2, list, lptr, lend, lp21)
Definition: tools_stripack.F90:1547
tools_stripack::bnodes
subroutine, public bnodes(n, list, lptr, lend, nodes, nb, na, nt)
Definition: tools_stripack.F90:421
tools_repro::small
logical function, public small(x, y)
Small value test.
Definition: tools_repro.F90:157
tools_stripack
STRIPACK routines.
Definition: tools_stripack.F90:11
tools_repro
Reproducibility functions.
Definition: tools_repro.F90:8
tools_stripack::jrand
subroutine jrand(n, ix, iy, iz, output)
Definition: tools_stripack.F90:1337
tools_repro::eq
logical function, public eq(x, y)
Equal test for reals.
Definition: tools_repro.F90:30
tools_stripack::trfind
subroutine, public trfind(nst, p, n, x, y, z, list, lptr, lend, b1, b2, b3, i1, i2, i3)
Definition: tools_stripack.F90:1765
tools_stripack::trmesh
subroutine, public trmesh(mpl, n, x, y, z, list, lptr, lend, lnew, near, next, dist, ier)
Definition: tools_stripack.F90:2579
tools_stripack::lstptr
subroutine lstptr(lpl, nb, list, lptr, output)
Definition: tools_stripack.F90:1466
tools_stripack::intadd
subroutine intadd(kk, i1, i2, i3, list, lptr, lend, lnew)
Definition: tools_stripack.F90:1127
tools_stripack::insert
subroutine insert(k, lp, list, lptr, lnew)
Definition: tools_stripack.F90:705
tools_stripack::addnod
subroutine, public addnod(mpl, nst, k, x, y, z, list, lptr, lend, lnew, ier)
Definition: tools_stripack.F90:25
tools_kinds
Kinds definition.
Definition: tools_kinds.F90:8
tools_stripack::det
subroutine det(x1, y1, z1, x2, y2, z2, x0, y0, z0, output)
Definition: tools_stripack.F90:670
type_mpl
MPI parameters derived type.
Definition: type_mpl.F90:8
tools_stripack::trlist
subroutine, public trlist(n, list, lptr, lend, nrow, nt, ltri, ier)
Definition: tools_stripack.F90:2279
type_mpl::mpl_type
Definition: type_mpl.F90:24
tools_stripack::left
subroutine left(x1, y1, z1, x2, y2, z2, x0, y0, z0, output)
Definition: tools_stripack.F90:1403
tools_kinds::kind_real
integer, parameter, public kind_real
Definition: tools_kinds.F90:18