SUBROUTINE ADDNOD (NST,K,X,Y,Z, LIST,LPTR,LEND, LNEW, IER) INTEGER NST, K, LIST(*), LPTR(*), LEND(K), LNEW, IER REAL X(K), Y(K), Z(K) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/08/99 ! ! This subroutine adds node K to a triangulation of the ! convex hull of nodes 1,...,K-1, producing a triangulation ! of the convex hull of nodes 1,...,K. ! ! The algorithm consists of the following steps: node K ! is located relative to the triangulation (TRFIND), its ! index is added to the data structure (INTADD or BDYADD), ! and a sequence of swaps (SWPTST and SWAP) are applied to ! the arcs opposite K so that all arcs incident on node K ! and opposite node K are locally optimal (satisfy the cir- ! cumcircle test). Thus, if a Delaunay triangulation is ! input, a Delaunay triangulation will result. ! ! ! On input: ! ! NST = Index of a node at which TRFIND begins its ! search. Search time depends on the proximity ! of this node to K. If NST < 1, the search is ! begun at node K-1. ! ! K = Nodal index (index for X, Y, Z, and LEND) of the ! new node to be added. K .GE. 4. ! ! X,Y,Z = Arrays of length .GE. K containing Car- ! tesian coordinates of the nodes. ! (X(I),Y(I),Z(I)) defines node I for ! I = 1,...,K. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LEND,LNEW = Data structure associated with ! the triangulation of nodes 1 ! to K-1. The array lengths are ! assumed to be large enough to ! add node K. Refer to Subrou- ! tine TRMESH. ! ! On output: ! ! LIST,LPTR,LEND,LNEW = Data structure updated with ! the addition of node K as the ! last entry unless IER .NE. 0 ! and IER .NE. -3, in which case ! the arrays are not altered. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = -1 if K is outside its valid range ! on input. ! IER = -2 if all nodes (including K) are col- ! linear (lie on a common geodesic). ! IER = L if nodes L and K coincide for some ! L < K. ! ! Modules required by ADDNOD: BDYADD, COVSPH, INSERT, ! INTADD, JRAND, LSTPTR, ! STORE, SWAP, SWPTST, ! TRFIND ! ! Intrinsic function called by ADDNOD: ABS ! !*********************************************************** ! INTEGER LSTPTR INTEGER I1, I2, I3, IO1, IO2, IN1, IST, KK, KM1, L, & & LP, LPF, LPO1, LPO1S LOGICAL SWPTST REAL B1, B2, B3, P(3) ! ! Local parameters: ! ! B1,B2,B3 = Unnormalized barycentric coordinates returned ! by TRFIND. ! I1,I2,I3 = Vertex indexes of a triangle containing K ! IN1 = Vertex opposite K: first neighbor of IO2 ! that precedes IO1. IN1,IO1,IO2 are in ! counterclockwise order. ! IO1,IO2 = Adjacent neighbors of K defining an arc to ! be tested for a swap ! IST = Index of node at which TRFIND begins its search ! KK = Local copy of K ! KM1 = K-1 ! L = Vertex index (I1, I2, or I3) returned in IER ! if node K coincides with a vertex ! LP = LIST pointer ! LPF = LIST pointer to the first neighbor of K ! LPO1 = LIST pointer to IO1 ! LPO1S = Saved value of LPO1 ! P = Cartesian coordinates of node K ! KK = K IF (KK .LT. 4) GO TO 3 ! ! Initialization: ! KM1 = KK - 1 IST = NST IF (IST .LT. 1) IST = KM1 P(1) = X(KK) P(2) = Y(KK) P(3) = Z(KK) ! ! Find a triangle (I1,I2,I3) containing K or the rightmost ! (I1) and leftmost (I2) visible boundary nodes as viewed ! from node K. ! CALL TRFIND (IST,P,KM1,X,Y,Z,LIST,LPTR,LEND, B1,B2,B3, I1,I2,I3) ! ! Test for collinear or duplicate nodes. ! IF (I1 .EQ. 0) GO TO 4 IF (I3 .NE. 0) THEN L = I1 IF (P(1) .EQ. X(L) .AND. P(2) .EQ. Y(L) .AND. & & P(3) .EQ. Z(L)) GO TO 5 L = I2 IF (P(1) .EQ. X(L) .AND. P(2) .EQ. Y(L) .AND. & & P(3) .EQ. Z(L)) GO TO 5 L = I3 IF (P(1) .EQ. X(L) .AND. P(2) .EQ. Y(L) .AND. & & P(3) .EQ. Z(L)) GO TO 5 CALL INTADD (KK,I1,I2,I3, LIST,LPTR,LEND,LNEW ) ELSE IF (I1 .NE. I2) THEN CALL BDYADD (KK,I1,I2, LIST,LPTR,LEND,LNEW ) ELSE CALL COVSPH (KK,I1, LIST,LPTR,LEND,LNEW ) ENDIF ENDIF IER = 0 ! ! Initialize variables for optimization of the ! triangulation. ! LP = LEND(KK) LPF = LPTR(LP) IO2 = LIST(LPF) LPO1 = LPTR(LPF) IO1 = ABS(LIST(LPO1)) ! ! Begin loop: find the node opposite K. ! 1 LP = LSTPTR(LEND(IO1),IO2,LIST,LPTR) IF (LIST(LP) .LT. 0) GO TO 2 LP = LPTR(LP) IN1 = ABS(LIST(LP)) ! ! Swap test: if a swap occurs, two new arcs are ! opposite K and must be tested. ! LPO1S = LPO1 IF ( .NOT. SWPTST(IN1,KK,IO1,IO2,X,Y,Z) ) GO TO 2 CALL SWAP (IN1,KK,IO1,IO2, LIST,LPTR,LEND, LPO1) IF (LPO1 .EQ. 0) THEN ! ! A swap is not possible because KK and IN1 are already ! adjacent. This error in SWPTST only occurs in the ! neutral case and when there are nearly duplicate ! nodes. ! LPO1 = LPO1S GO TO 2 ENDIF IO1 = IN1 GO TO 1 ! ! No swap occurred. Test for termination and reset ! IO2 and IO1. ! 2 IF (LPO1 .EQ. LPF .OR. LIST(LPO1) .LT. 0) RETURN IO2 = IO1 LPO1 = LPTR(LPO1) IO1 = ABS(LIST(LPO1)) GO TO 1 ! ! KK < 4. ! 3 IER = -1 RETURN ! ! All nodes are collinear. ! 4 IER = -2 RETURN ! ! Nodes L and K coincide. ! 5 IER = L RETURN END REAL FUNCTION AREAS (V1,V2,V3) REAL V1(3), V2(3), V3(3) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 09/18/90 ! ! This function returns the area of a spherical triangle ! on the unit sphere. ! ! ! On input: ! ! V1,V2,V3 = Arrays of length 3 containing the Carte- ! sian coordinates of unit vectors (the ! three triangle vertices in any order). ! These vectors, if nonzero, are implicitly ! scaled to have length 1. ! ! Input parameters are not altered by this function. ! ! On output: ! ! AREAS = Area of the spherical triangle defined by ! V1, V2, and V3 in the range 0 to 2*PI (the ! area of a hemisphere). AREAS = 0 (or 2*PI) ! if and only if V1, V2, and V3 lie in (or ! close to) a plane containing the origin. ! ! Modules required by AREAS: None ! ! Intrinsic functions called by AREAS: ACOS, DBLE, REAL, ! SQRT ! !*********************************************************** ! DOUBLE PRECISION A1, A2, A3, CA1, CA2, CA3, DV1(3), & & DV2(3), DV3(3), S12, S23, S31, & & U12(3), U23(3), U31(3) INTEGER I ! ! Local parameters: ! ! A1,A2,A3 = Interior angles of the spherical triangle ! CA1,CA2,CA3 = cos(A1), cos(A2), and cos(A3), respectively ! DV1,DV2,DV3 = Double Precision copies of V1, V2, and V3 ! I = DO-loop index and index for Uij ! S12,S23,S31 = Sum of squared components of U12, U23, U31 ! U12,U23,U31 = Unit normal vectors to the planes defined by ! pairs of triangle vertices ! DO 1 I = 1,3 DV1(I) = DBLE(V1(I)) DV2(I) = DBLE(V2(I)) DV3(I) = DBLE(V3(I)) 1 CONTINUE ! ! Compute cross products Uij = Vi X Vj. ! U12(1) = DV1(2)*DV2(3) - DV1(3)*DV2(2) U12(2) = DV1(3)*DV2(1) - DV1(1)*DV2(3) U12(3) = DV1(1)*DV2(2) - DV1(2)*DV2(1) ! U23(1) = DV2(2)*DV3(3) - DV2(3)*DV3(2) U23(2) = DV2(3)*DV3(1) - DV2(1)*DV3(3) U23(3) = DV2(1)*DV3(2) - DV2(2)*DV3(1) ! U31(1) = DV3(2)*DV1(3) - DV3(3)*DV1(2) U31(2) = DV3(3)*DV1(1) - DV3(1)*DV1(3) U31(3) = DV3(1)*DV1(2) - DV3(2)*DV1(1) ! ! Normalize Uij to unit vectors. ! S12 = 0.D0 S23 = 0.D0 S31 = 0.D0 DO 2 I = 1,3 S12 = S12 + U12(I)*U12(I) S23 = S23 + U23(I)*U23(I) S31 = S31 + U31(I)*U31(I) 2 CONTINUE ! ! Test for a degenerate triangle associated with collinear ! vertices. ! IF (S12 .EQ. 0.D0 .OR. S23 .EQ. 0.D0 .OR. & & S31 .EQ. 0.D0) THEN AREAS = 0. RETURN ENDIF S12 = SQRT(S12) S23 = SQRT(S23) S31 = SQRT(S31) DO 3 I = 1,3 U12(I) = U12(I)/S12 U23(I) = U23(I)/S23 U31(I) = U31(I)/S31 3 CONTINUE ! ! Compute interior angles Ai as the dihedral angles between ! planes: ! CA1 = cos(A1) = - ! CA2 = cos(A2) = - ! CA3 = cos(A3) = - ! CA1 = -U12(1)*U31(1)-U12(2)*U31(2)-U12(3)*U31(3) CA2 = -U23(1)*U12(1)-U23(2)*U12(2)-U23(3)*U12(3) CA3 = -U31(1)*U23(1)-U31(2)*U23(2)-U31(3)*U23(3) IF (CA1 .LT. -1.D0) CA1 = -1.D0 IF (CA1 .GT. 1.D0) CA1 = 1.D0 IF (CA2 .LT. -1.D0) CA2 = -1.D0 IF (CA2 .GT. 1.D0) CA2 = 1.D0 IF (CA3 .LT. -1.D0) CA3 = -1.D0 IF (CA3 .GT. 1.D0) CA3 = 1.D0 A1 = ACOS(CA1) A2 = ACOS(CA2) A3 = ACOS(CA3) ! ! Compute AREAS = A1 + A2 + A3 - PI. ! AREAS = REAL(A1 + A2 + A3 - ACOS(-1.D0)) IF (AREAS .LT. 0.) AREAS = 0. RETURN END SUBROUTINE BDYADD (KK,I1,I2, LIST,LPTR,LEND,LNEW ) INTEGER KK, I1, I2, LIST(*), LPTR(*), LEND(*), LNEW ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/11/96 ! ! This subroutine adds a boundary node to a triangulation ! of a set of KK-1 points on the unit sphere. The data ! structure is updated with the insertion of node KK, but no ! optimization is performed. ! ! This routine is identical to the similarly named routine ! in TRIPACK. ! ! ! On input: ! ! KK = Index of a node to be connected to the sequence ! of all visible boundary nodes. KK .GE. 1 and ! KK must not be equal to I1 or I2. ! ! I1 = First (rightmost as viewed from KK) boundary ! node in the triangulation that is visible from ! node KK (the line segment KK-I1 intersects no ! arcs. ! ! I2 = Last (leftmost) boundary node that is visible ! from node KK. I1 and I2 may be determined by ! Subroutine TRFIND. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LEND,LNEW = Triangulation data structure ! created by Subroutine TRMESH. ! Nodes I1 and I2 must be in- ! cluded in the triangulation. ! ! On output: ! ! LIST,LPTR,LEND,LNEW = Data structure updated with ! the addition of node KK. Node ! KK is connected to I1, I2, and ! all boundary nodes in between. ! ! Module required by BDYADD: INSERT ! !*********************************************************** ! INTEGER K, LP, LSAV, N1, N2, NEXT, NSAV ! ! Local parameters: ! ! K = Local copy of KK ! LP = LIST pointer ! LSAV = LIST pointer ! N1,N2 = Local copies of I1 and I2, respectively ! NEXT = Boundary node visible from K ! NSAV = Boundary node visible from K ! K = KK N1 = I1 N2 = I2 ! ! Add K as the last neighbor of N1. ! LP = LEND(N1) LSAV = LPTR(LP) LPTR(LP) = LNEW LIST(LNEW) = -K LPTR(LNEW) = LSAV LEND(N1) = LNEW LNEW = LNEW + 1 NEXT = -LIST(LP) LIST(LP) = NEXT NSAV = NEXT ! ! Loop on the remaining boundary nodes between N1 and N2, ! adding K as the first neighbor. ! 1 LP = LEND(NEXT) CALL INSERT (K,LP, LIST,LPTR,LNEW ) IF (NEXT .EQ. N2) GO TO 2 NEXT = -LIST(LP) LIST(LP) = NEXT GO TO 1 ! ! Add the boundary nodes between N1 and N2 as neighbors ! of node K. ! 2 LSAV = LNEW LIST(LNEW) = N1 LPTR(LNEW) = LNEW + 1 LNEW = LNEW + 1 NEXT = NSAV ! 3 IF (NEXT .EQ. N2) GO TO 4 LIST(LNEW) = NEXT LPTR(LNEW) = LNEW + 1 LNEW = LNEW + 1 LP = LEND(NEXT) NEXT = LIST(LP) GO TO 3 ! 4 LIST(LNEW) = -N2 LPTR(LNEW) = LSAV LEND(K) = LNEW LNEW = LNEW + 1 RETURN END SUBROUTINE BNODES (N,LIST,LPTR,LEND, NODES,NB,NA,NT) INTEGER N, LIST(*), LPTR(*), LEND(N), NODES(*), NB, NA, NT ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 06/26/96 ! ! Given a triangulation of N nodes on the unit sphere ! created by Subroutine TRMESH, this subroutine returns an ! array containing the indexes (if any) of the counterclock- ! wise-ordered sequence of boundary nodes -- the nodes on ! the boundary of the convex hull of the set of nodes. (The ! boundary is empty if the nodes do not lie in a single ! hemisphere.) The numbers of boundary nodes, arcs, and ! triangles are also returned. ! ! ! On input: ! ! N = Number of nodes in the triangulation. N .GE. 3. ! ! LIST,LPTR,LEND = Data structure defining the trian- ! gulation. Refer to Subroutine ! TRMESH. ! ! The above parameters are not altered by this routine. ! ! NODES = Integer array of length at least NB ! (NB .LE. N). ! ! On output: ! ! NODES = Ordered sequence of boundary node indexes ! in the range 1 to N (in the first NB loca- ! tions). ! ! NB = Number of boundary nodes. ! ! NA,NT = Number of arcs and triangles, respectively, ! in the triangulation. ! ! Modules required by BNODES: None ! !*********************************************************** ! INTEGER K, LP, N0, NN, NST ! ! Local parameters: ! ! K = NODES index ! LP = LIST pointer ! N0 = Boundary node to be added to NODES ! NN = Local copy of N ! NST = First element of nodes (arbitrarily chosen to be ! the one with smallest index) ! NN = N ! ! Search for a boundary node. ! DO 1 NST = 1,NN LP = LEND(NST) IF (LIST(LP) .LT. 0) GO TO 2 1 CONTINUE ! ! The triangulation contains no boundary nodes. ! NB = 0 NA = 3*(NN-2) NT = 2*(NN-2) RETURN ! ! NST is the first boundary node encountered. Initialize ! for traversal of the boundary. ! 2 NODES(1) = NST K = 1 N0 = NST ! ! Traverse the boundary in counterclockwise order. ! 3 LP = LEND(N0) LP = LPTR(LP) N0 = LIST(LP) IF (N0 .EQ. NST) GO TO 4 K = K + 1 NODES(K) = N0 GO TO 3 ! ! Store the counts. ! 4 NB = K NT = 2*N - NB - 2 NA = NT + N - 1 RETURN END SUBROUTINE CIRCUM (V1,V2,V3, C,IER) INTEGER IER REAL V1(3), V2(3), V3(3), C(3) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 06/29/95 ! ! This subroutine returns the circumcenter of a spherical ! triangle on the unit sphere: the point on the sphere sur- ! face that is equally distant from the three triangle ! vertices and lies in the same hemisphere, where distance ! is taken to be arc-length on the sphere surface. ! ! ! On input: ! ! V1,V2,V3 = Arrays of length 3 containing the Carte- ! sian coordinates of the three triangle ! vertices (unit vectors) in CCW order. ! ! The above parameters are not altered by this routine. ! ! C = Array of length 3. ! ! On output: ! ! C = Cartesian coordinates of the circumcenter unless ! IER > 0, in which case C is not defined. C = ! (V2-V1) X (V3-V1) normalized to a unit vector. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if V1, V2, and V3 lie on a common ! line: (V2-V1) X (V3-V1) = 0. ! (The vertices are not tested for validity.) ! ! Modules required by CIRCUM: None ! ! Intrinsic function called by CIRCUM: SQRT ! !*********************************************************** ! INTEGER I REAL CNORM, CU(3), E1(3), E2(3) ! ! Local parameters: ! ! CNORM = Norm of CU: used to compute C ! CU = Scalar multiple of C: E1 X E2 ! E1,E2 = Edges of the underlying planar triangle: ! V2-V1 and V3-V1, respectively ! I = DO-loop index ! DO 1 I = 1,3 E1(I) = V2(I) - V1(I) E2(I) = V3(I) - V1(I) 1 CONTINUE ! ! Compute CU = E1 X E2 and CNORM**2. ! CU(1) = E1(2)*E2(3) - E1(3)*E2(2) CU(2) = E1(3)*E2(1) - E1(1)*E2(3) CU(3) = E1(1)*E2(2) - E1(2)*E2(1) CNORM = CU(1)*CU(1) + CU(2)*CU(2) + CU(3)*CU(3) ! ! The vertices lie on a common line if and only if CU is ! the zero vector. ! IF (CNORM .NE. 0.) THEN ! ! No error: compute C. ! CNORM = SQRT(CNORM) DO 2 I = 1,3 C(I) = CU(I)/CNORM 2 CONTINUE IER = 0 ELSE ! ! CU = 0. ! IER = 1 ENDIF RETURN END SUBROUTINE COVSPH (KK,N0, LIST,LPTR,LEND,LNEW ) INTEGER KK, N0, LIST(*), LPTR(*), LEND(*), LNEW ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/17/96 ! ! This subroutine connects an exterior node KK to all ! boundary nodes of a triangulation of KK-1 points on the ! unit sphere, producing a triangulation that covers the ! sphere. The data structure is updated with the addition ! of node KK, but no optimization is performed. All boun- ! dary nodes must be visible from node KK. ! ! ! On input: ! ! KK = Index of the node to be connected to the set of ! all boundary nodes. KK .GE. 4. ! ! N0 = Index of a boundary node (in the range 1 to ! KK-1). N0 may be determined by Subroutine ! TRFIND. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LEND,LNEW = Triangulation data structure ! created by Subroutine TRMESH. ! Node N0 must be included in ! the triangulation. ! ! On output: ! ! LIST,LPTR,LEND,LNEW = Data structure updated with ! the addition of node KK as the ! last entry. The updated ! triangulation contains no ! boundary nodes. ! ! Module required by COVSPH: INSERT ! !*********************************************************** ! INTEGER K, LP, LSAV, NEXT, NST ! ! Local parameters: ! ! K = Local copy of KK ! LP = LIST pointer ! LSAV = LIST pointer ! NEXT = Boundary node visible from K ! NST = Local copy of N0 ! K = KK NST = N0 ! ! Traverse the boundary in clockwise order, inserting K as ! the first neighbor of each boundary node, and converting ! the boundary node to an interior node. ! NEXT = NST 1 LP = LEND(NEXT) CALL INSERT (K,LP, LIST,LPTR,LNEW ) NEXT = -LIST(LP) LIST(LP) = NEXT IF (NEXT .NE. NST) GO TO 1 ! ! Traverse the boundary again, adding each node to K's ! adjacency list. ! LSAV = LNEW 2 LP = LEND(NEXT) LIST(LNEW) = NEXT LPTR(LNEW) = LNEW + 1 LNEW = LNEW + 1 NEXT = LIST(LP) IF (NEXT .NE. NST) GO TO 2 ! LPTR(LNEW-1) = LSAV LEND(K) = LNEW - 1 RETURN END SUBROUTINE CRLIST (N,NCOL,X,Y,Z,LIST,LEND, LPTR,LNEW, & & LTRI, LISTC,NB,XC,YC,ZC,RC,IER) INTEGER N, NCOL, LIST(*), LEND(N), LPTR(*), LNEW, & & LTRI(6,NCOL), LISTC(*), NB, IER REAL X(N), Y(N), Z(N), XC(*), YC(*), ZC(*), RC(*) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/05/98 ! ! Given a Delaunay triangulation of nodes on the surface ! of the unit sphere, this subroutine returns the set of ! triangle circumcenters corresponding to Voronoi vertices, ! along with the circumradii and a list of triangle indexes ! LISTC stored in one-to-one correspondence with LIST/LPTR ! entries. ! ! A triangle circumcenter is the point (unit vector) lying ! at the same angular distance from the three vertices and ! contained in the same hemisphere as the vertices. (Note ! that the negative of a circumcenter is also equidistant ! from the vertices.) If the triangulation covers the sur- ! face, the Voronoi vertices are the circumcenters of the ! triangles in the Delaunay triangulation. LPTR, LEND, and ! LNEW are not altered in this case. ! ! On the other hand, if the nodes are contained in a sin- ! gle hemisphere, the triangulation is implicitly extended ! to the entire surface by adding pseudo-arcs (of length ! greater than 180 degrees) between boundary nodes forming ! pseudo-triangles whose 'circumcenters' are included in the ! list. This extension to the triangulation actually con- ! sists of a triangulation of the set of boundary nodes in ! which the swap test is reversed (a non-empty circumcircle ! test). The negative circumcenters are stored as the ! pseudo-triangle 'circumcenters'. LISTC, LPTR, LEND, and ! LNEW contain a data structure corresponding to the ex- ! tended triangulation (Voronoi diagram), but LIST is not ! altered in this case. Thus, if it is necessary to retain ! the original (unextended) triangulation data structure, ! copies of LPTR and LNEW must be saved before calling this ! routine. ! ! ! On input: ! ! N = Number of nodes in the triangulation. N .GE. 3. ! Note that, if N = 3, there are only two Voronoi ! vertices separated by 180 degrees, and the ! Voronoi regions are not well defined. ! ! NCOL = Number of columns reserved for LTRI. This ! must be at least NB-2, where NB is the number ! of boundary nodes. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes (unit vectors). ! ! LIST = Integer array containing the set of adjacency ! lists. Refer to Subroutine TRMESH. ! ! LEND = Set of pointers to ends of adjacency lists. ! Refer to Subroutine TRMESH. ! ! The above parameters are not altered by this routine. ! ! LPTR = Array of pointers associated with LIST. Re- ! fer to Subroutine TRMESH. ! ! LNEW = Pointer to the first empty location in LIST ! and LPTR (list length plus one). ! ! LTRI = Integer work space array dimensioned 6 by ! NCOL, or unused dummy parameter if NB = 0. ! ! LISTC = Integer array of length at least 3*NT, where ! NT = 2*N-4 is the number of triangles in the ! triangulation (after extending it to cover ! the entire surface if necessary). ! ! XC,YC,ZC,RC = Arrays of length NT = 2*N-4. ! ! On output: ! ! LPTR = Array of pointers associated with LISTC: ! updated for the addition of pseudo-triangles ! if the original triangulation contains ! boundary nodes (NB > 0). ! ! LNEW = Pointer to the first empty location in LISTC ! and LPTR (list length plus one). LNEW is not ! altered if NB = 0. ! ! LTRI = Triangle list whose first NB-2 columns con- ! tain the indexes of a clockwise-ordered ! sequence of vertices (first three rows) ! followed by the LTRI column indexes of the ! triangles opposite the vertices (or 0 ! denoting the exterior region) in the last ! three rows. This array is not generally of ! any use. ! ! LISTC = Array containing triangle indexes (indexes ! to XC, YC, ZC, and RC) stored in 1-1 corres- ! pondence with LIST/LPTR entries (or entries ! that would be stored in LIST for the ! extended triangulation): the index of tri- ! angle (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). ! ! NB = Number of boundary nodes unless IER = 1. ! ! XC,YC,ZC = Arrays containing the Cartesian coordi- ! nates of the triangle circumcenters ! (Voronoi vertices). XC(I)**2 + YC(I)**2 ! + ZC(I)**2 = 1. The first NB-2 entries ! correspond to pseudo-triangles if NB > 0. ! ! RC = Array containing circumradii (the arc lengths ! or angles between the circumcenters and associ- ! ated triangle vertices) in 1-1 correspondence ! with circumcenters. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if N < 3. ! IER = 2 if NCOL < NB-2. ! IER = 3 if a triangle is degenerate (has ver- ! tices lying on a common geodesic). ! ! Modules required by CRLIST: CIRCUM, LSTPTR, SWPTST ! ! Intrinsic functions called by CRLIST: ABS, ACOS ! !*********************************************************** ! INTEGER LSTPTR INTEGER I1, I2, I3, I4, IERR, KT, KT1, KT2, KT11, & & KT12, KT21, KT22, LP, LPL, LPN, N0, N1, N2, & & N3, N4, NM2, NN, NT LOGICAL SWPTST LOGICAL SWP REAL C(3), T, V1(3), V2(3), V3(3) ! ! Local parameters: ! ! C = Circumcenter returned by Subroutine CIRCUM ! I1,I2,I3 = Permutation of (1,2,3): LTRI row indexes ! I4 = LTRI row index in the range 1 to 3 ! IERR = Error flag for calls to CIRCUM ! KT = Triangle index ! KT1,KT2 = Indexes of a pair of adjacent pseudo-triangles ! KT11,KT12 = Indexes of the pseudo-triangles opposite N1 ! and N2 as vertices of KT1 ! KT21,KT22 = Indexes of the pseudo-triangles opposite N1 ! and N2 as vertices of KT2 ! LP,LPN = LIST pointers ! LPL = LIST pointer of the last neighbor of N1 ! N0 = Index of the first boundary node (initial ! value of N1) in the loop on boundary nodes ! used to store the pseudo-triangle indexes ! in LISTC ! N1,N2,N3 = Nodal indexes defining a triangle (CCW order) ! or pseudo-triangle (clockwise order) ! N4 = Index of the node opposite N2 -> N1 ! NM2 = N-2 ! NN = Local copy of N ! NT = Number of pseudo-triangles: NB-2 ! SWP = Logical variable set to TRUE in each optimiza- ! tion loop (loop on pseudo-arcs) iff a swap ! is performed ! V1,V2,V3 = Vertices of triangle KT = (N1,N2,N3) sent to ! Subroutine CIRCUM ! NN = N NB = 0 NT = 0 IF (NN .LT. 3) GO TO 21 ! ! Search for a boundary node N1. ! DO 1 N1 = 1,NN LP = LEND(N1) IF (LIST(LP) .LT. 0) GO TO 2 1 CONTINUE ! ! The triangulation already covers the sphere. ! GO TO 9 ! ! There are NB .GE. 3 boundary nodes. Add NB-2 pseudo- ! triangles (N1,N2,N3) by connecting N3 to the NB-3 ! boundary nodes to which it is not already adjacent. ! ! Set N3 and N2 to the first and last neighbors, ! respectively, of N1. ! 2 N2 = -LIST(LP) LP = LPTR(LP) N3 = LIST(LP) ! ! Loop on boundary arcs N1 -> N2 in clockwise order, ! storing triangles (N1,N2,N3) in column NT of LTRI ! along with the indexes of the triangles opposite ! the vertices. ! 3 NT = NT + 1 IF (NT .LE. NCOL) THEN LTRI(1,NT) = N1 LTRI(2,NT) = N2 LTRI(3,NT) = N3 LTRI(4,NT) = NT + 1 LTRI(5,NT) = NT - 1 LTRI(6,NT) = 0 ENDIF N1 = N2 LP = LEND(N1) N2 = -LIST(LP) IF (N2 .NE. N3) GO TO 3 ! NB = NT + 2 IF (NCOL .LT. NT) GO TO 22 LTRI(4,NT) = 0 IF (NT .EQ. 1) GO TO 7 ! ! Optimize the exterior triangulation (set of pseudo- ! triangles) by applying swaps to the pseudo-arcs N1-N2 ! (pairs of adjacent pseudo-triangles KT1 and KT2 > KT1). ! The loop on pseudo-arcs is repeated until no swaps are ! performed. ! 4 SWP = .FALSE. DO 6 KT1 = 1,NT-1 DO 5 I3 = 1,3 KT2 = LTRI(I3+3,KT1) IF (KT2 .LE. KT1) GO TO 5 ! ! The LTRI row indexes (I1,I2,I3) of triangle KT1 = ! (N1,N2,N3) are a cyclical permutation of (1,2,3). ! IF (I3 .EQ. 1) THEN I1 = 2 I2 = 3 ELSEIF (I3 .EQ. 2) THEN I1 = 3 I2 = 1 ELSE I1 = 1 I2 = 2 ENDIF N1 = LTRI(I1,KT1) N2 = LTRI(I2,KT1) N3 = LTRI(I3,KT1) ! ! KT2 = (N2,N1,N4) for N4 = LTRI(I,KT2), where ! LTRI(I+3,KT2) = KT1. ! IF (LTRI(4,KT2) .EQ. KT1) THEN I4 = 1 ELSEIF (LTRI(5,KT2) .EQ. KT1) THEN I4 = 2 ELSE I4 = 3 ENDIF N4 = LTRI(I4,KT2) ! ! The empty circumcircle test is reversed for the pseudo- ! triangles. The reversal is implicit in the clockwise ! ordering of the vertices. ! IF ( .NOT. SWPTST(N1,N2,N3,N4,X,Y,Z) ) GO TO 5 ! ! Swap arc N1-N2 for N3-N4. KTij is the triangle opposite ! Nj as a vertex of KTi. ! SWP = .TRUE. KT11 = LTRI(I1+3,KT1) KT12 = LTRI(I2+3,KT1) IF (I4 .EQ. 1) THEN I2 = 2 I1 = 3 ELSEIF (I4 .EQ. 2) THEN I2 = 3 I1 = 1 ELSE I2 = 1 I1 = 2 ENDIF KT21 = LTRI(I1+3,KT2) KT22 = LTRI(I2+3,KT2) LTRI(1,KT1) = N4 LTRI(2,KT1) = N3 LTRI(3,KT1) = N1 LTRI(4,KT1) = KT12 LTRI(5,KT1) = KT22 LTRI(6,KT1) = KT2 LTRI(1,KT2) = N3 LTRI(2,KT2) = N4 LTRI(3,KT2) = N2 LTRI(4,KT2) = KT21 LTRI(5,KT2) = KT11 LTRI(6,KT2) = KT1 ! ! Correct the KT11 and KT22 entries that changed. ! IF (KT11 .NE. 0) THEN I4 = 4 IF (LTRI(4,KT11) .NE. KT1) THEN I4 = 5 IF (LTRI(5,KT11) .NE. KT1) I4 = 6 ENDIF LTRI(I4,KT11) = KT2 ENDIF IF (KT22 .NE. 0) THEN I4 = 4 IF (LTRI(4,KT22) .NE. KT2) THEN I4 = 5 IF (LTRI(5,KT22) .NE. KT2) I4 = 6 ENDIF LTRI(I4,KT22) = KT1 ENDIF 5 CONTINUE 6 CONTINUE IF (SWP) GO TO 4 ! ! Compute and store the negative circumcenters and radii of ! the pseudo-triangles in the first NT positions. ! 7 DO 8 KT = 1,NT N1 = LTRI(1,KT) N2 = LTRI(2,KT) N3 = LTRI(3,KT) V1(1) = X(N1) V1(2) = Y(N1) V1(3) = Z(N1) V2(1) = X(N2) V2(2) = Y(N2) V2(3) = Z(N2) V3(1) = X(N3) V3(2) = Y(N3) V3(3) = Z(N3) CALL CIRCUM (V1,V2,V3, C,IERR) IF (IERR .NE. 0) GO TO 23 ! ! Store the negative circumcenter and radius (computed ! from ). ! XC(KT) = C(1) YC(KT) = C(2) ZC(KT) = C(3) T = V1(1)*C(1) + V1(2)*C(2) + V1(3)*C(3) IF (T .LT. -1.0) T = -1.0 IF (T .GT. 1.0) T = 1.0 RC(KT) = ACOS(T) 8 CONTINUE ! ! Compute and store the circumcenters and radii of the ! actual triangles in positions KT = NT+1, NT+2, ... ! Also, store the triangle indexes KT in the appropriate ! LISTC positions. ! 9 KT = NT ! ! Loop on nodes N1. ! NM2 = NN - 2 DO 12 N1 = 1,NM2 LPL = LEND(N1) LP = LPL N3 = LIST(LP) ! ! Loop on adjacent neighbors N2,N3 of N1 for which N2 > N1 ! and N3 > N1. ! 10 LP = LPTR(LP) N2 = N3 N3 = ABS(LIST(LP)) IF (N2 .LE. N1 .OR. N3 .LE. N1) GO TO 11 KT = KT + 1 ! ! Compute the circumcenter C of triangle KT = (N1,N2,N3). ! V1(1) = X(N1) V1(2) = Y(N1) V1(3) = Z(N1) V2(1) = X(N2) V2(2) = Y(N2) V2(3) = Z(N2) V3(1) = X(N3) V3(2) = Y(N3) V3(3) = Z(N3) CALL CIRCUM (V1,V2,V3, C,IERR) IF (IERR .NE. 0) GO TO 23 ! ! Store the circumcenter, radius and triangle index. ! XC(KT) = C(1) YC(KT) = C(2) ZC(KT) = C(3) T = V1(1)*C(1) + V1(2)*C(2) + V1(3)*C(3) IF (T .LT. -1.0) T = -1.0 IF (T .GT. 1.0) T = 1.0 RC(KT) = ACOS(T) ! ! Store KT in LISTC(LPN), where Abs(LIST(LPN)) is the ! index of N2 as a neighbor of N1, N3 as a neighbor ! of N2, and N1 as a neighbor of N3. ! LPN = LSTPTR(LPL,N2,LIST,LPTR) LISTC(LPN) = KT LPN = LSTPTR(LEND(N2),N3,LIST,LPTR) LISTC(LPN) = KT LPN = LSTPTR(LEND(N3),N1,LIST,LPTR) LISTC(LPN) = KT 11 IF (LP .NE. LPL) GO TO 10 12 CONTINUE IF (NT .EQ. 0) GO TO 20 ! ! Store the first NT triangle indexes in LISTC. ! ! Find a boundary triangle KT1 = (N1,N2,N3) with a ! boundary arc opposite N3. ! KT1 = 0 13 KT1 = KT1 + 1 IF (LTRI(4,KT1) .EQ. 0) THEN I1 = 2 I2 = 3 I3 = 1 GO TO 14 ELSEIF (LTRI(5,KT1) .EQ. 0) THEN I1 = 3 I2 = 1 I3 = 2 GO TO 14 ELSEIF (LTRI(6,KT1) .EQ. 0) THEN I1 = 1 I2 = 2 I3 = 3 GO TO 14 ENDIF GO TO 13 14 N1 = LTRI(I1,KT1) N0 = N1 ! ! Loop on boundary nodes N1 in CCW order, storing the ! indexes of the clockwise-ordered sequence of triangles ! that contain N1. The first triangle overwrites the ! last neighbor position, and the remaining triangles, ! if any, are appended to N1's adjacency list. ! ! A pointer to the first neighbor of N1 is saved in LPN. ! 15 LP = LEND(N1) LPN = LPTR(LP) LISTC(LP) = KT1 ! ! Loop on triangles KT2 containing N1. ! 16 KT2 = LTRI(I2+3,KT1) IF (KT2 .NE. 0) THEN ! ! Append KT2 to N1's triangle list. ! LPTR(LP) = LNEW LP = LNEW LISTC(LP) = KT2 LNEW = LNEW + 1 ! ! Set KT1 to KT2 and update (I1,I2,I3) such that ! LTRI(I1,KT1) = N1. ! KT1 = KT2 IF (LTRI(1,KT1) .EQ. N1) THEN I1 = 1 I2 = 2 I3 = 3 ELSEIF (LTRI(2,KT1) .EQ. N1) THEN I1 = 2 I2 = 3 I3 = 1 ELSE I1 = 3 I2 = 1 I3 = 2 ENDIF GO TO 16 ENDIF ! ! Store the saved first-triangle pointer in LPTR(LP), set ! N1 to the next boundary node, test for termination, ! and permute the indexes: the last triangle containing ! a boundary node is the first triangle containing the ! next boundary node. ! LPTR(LP) = LPN N1 = LTRI(I3,KT1) IF (N1 .NE. N0) THEN I4 = I3 I3 = I2 I2 = I1 I1 = I4 GO TO 15 ENDIF ! ! No errors encountered. ! 20 IER = 0 RETURN ! ! N < 3. ! 21 IER = 1 RETURN ! ! Insufficient space reserved for LTRI. ! 22 IER = 2 RETURN ! ! Error flag returned by CIRCUM: KT indexes a null triangle. ! 23 IER = 3 RETURN END SUBROUTINE DELARC (N,IO1,IO2, LIST,LPTR,LEND, LNEW, IER) INTEGER N, IO1, IO2, LIST(*), LPTR(*), LEND(N), LNEW, IER ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/17/96 ! ! This subroutine deletes a boundary arc from a triangula- ! tion. It may be used to remove a null triangle from the ! convex hull boundary. Note, however, that if the union of ! triangles is rendered nonconvex, Subroutines DELNOD, EDGE, ! and TRFIND (and hence ADDNOD) may fail. Also, Function ! NEARND should not be called following an arc deletion. ! ! This routine is identical to the similarly named routine ! in TRIPACK. ! ! ! On input: ! ! N = Number of nodes in the triangulation. N .GE. 4. ! ! IO1,IO2 = Indexes (in the range 1 to N) of a pair of ! adjacent boundary nodes defining the arc ! to be removed. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LEND,LNEW = Triangulation data structure ! created by Subroutine TRMESH. ! ! On output: ! ! LIST,LPTR,LEND,LNEW = Data structure updated with ! the removal of arc IO1-IO2 ! unless IER > 0. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if N, IO1, or IO2 is outside its valid ! range, or IO1 = IO2. ! IER = 2 if IO1-IO2 is not a boundary arc. ! IER = 3 if the node opposite IO1-IO2 is al- ! ready a boundary node, and thus IO1 ! or IO2 has only two neighbors or a ! deletion would result in two triangu- ! lations sharing a single node. ! IER = 4 if one of the nodes is a neighbor of ! the other, but not vice versa, imply- ! ing an invalid triangulation data ! structure. ! ! Module required by DELARC: DELNB, LSTPTR ! ! Intrinsic function called by DELARC: ABS ! !*********************************************************** ! INTEGER LSTPTR INTEGER LP, LPH, LPL, N1, N2, N3 ! ! Local parameters: ! ! LP = LIST pointer ! LPH = LIST pointer or flag returned by DELNB ! LPL = Pointer to the last neighbor of N1, N2, or N3 ! N1,N2,N3 = Nodal indexes of a triangle such that N1->N2 ! is the directed boundary edge associated ! with IO1-IO2 ! N1 = IO1 N2 = IO2 ! ! Test for errors, and set N1->N2 to the directed boundary ! edge associated with IO1-IO2: (N1,N2,N3) is a triangle ! for some N3. ! IF (N .LT. 4 .OR. N1 .LT. 1 .OR. N1 .GT. N .OR. & & N2 .LT. 1 .OR. N2 .GT. N .OR. N1 .EQ. N2) THEN IER = 1 RETURN ENDIF ! LPL = LEND(N2) IF (-LIST(LPL) .NE. N1) THEN N1 = N2 N2 = IO1 LPL = LEND(N2) IF (-LIST(LPL) .NE. N1) THEN IER = 2 RETURN ENDIF ENDIF ! ! Set N3 to the node opposite N1->N2 (the second neighbor ! of N1), and test for error 3 (N3 already a boundary ! node). ! LPL = LEND(N1) LP = LPTR(LPL) LP = LPTR(LP) N3 = ABS(LIST(LP)) LPL = LEND(N3) IF (LIST(LPL) .LE. 0) THEN IER = 3 RETURN ENDIF ! ! Delete N2 as a neighbor of N1, making N3 the first ! neighbor, and test for error 4 (N2 not a neighbor ! of N1). Note that previously computed pointers may ! no longer be valid following the call to DELNB. ! CALL DELNB (N1,N2,N, LIST,LPTR,LEND,LNEW, LPH) IF (LPH .LT. 0) THEN IER = 4 RETURN ENDIF ! ! Delete N1 as a neighbor of N2, making N3 the new last ! neighbor. ! CALL DELNB (N2,N1,N, LIST,LPTR,LEND,LNEW, LPH) ! ! Make N3 a boundary node with first neighbor N2 and last ! neighbor N1. ! LP = LSTPTR(LEND(N3),N1,LIST,LPTR) LEND(N3) = LP LIST(LP) = -N1 ! ! No errors encountered. ! IER = 0 RETURN END SUBROUTINE DELNB (N0,NB,N, LIST,LPTR,LEND,LNEW, LPH) INTEGER N0, NB, N, LIST(*), LPTR(*), LEND(N), LNEW, LPH ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/29/98 ! ! This subroutine deletes a neighbor NB from the adjacency ! list of node N0 (but N0 is not deleted from the adjacency ! list of NB) and, if NB is a boundary node, makes N0 a ! boundary node. For pointer (LIST index) LPH to NB as a ! neighbor of N0, the empty LIST,LPTR location LPH is filled ! in with the values at LNEW-1, pointer LNEW-1 (in LPTR and ! possibly in LEND) is changed to LPH, and LNEW is decremen- ! ted. This requires a search of LEND and LPTR entailing an ! expected operation count of O(N). ! ! This routine is identical to the similarly named routine ! in TRIPACK. ! ! ! On input: ! ! N0,NB = Indexes, in the range 1 to N, of a pair of ! nodes such that NB is a neighbor of N0. ! (N0 need not be a neighbor of NB.) ! ! N = Number of nodes in the triangulation. N .GE. 3. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LEND,LNEW = Data structure defining the ! triangulation. ! ! On output: ! ! LIST,LPTR,LEND,LNEW = Data structure updated with ! the removal of NB from the ad- ! jacency list of N0 unless ! LPH < 0. ! ! LPH = List pointer to the hole (NB as a neighbor of ! N0) filled in by the values at LNEW-1 or error ! indicator: ! LPH > 0 if no errors were encountered. ! LPH = -1 if N0, NB, or N is outside its valid ! range. ! LPH = -2 if NB is not a neighbor of N0. ! ! Modules required by DELNB: None ! ! Intrinsic function called by DELNB: ABS ! !*********************************************************** ! INTEGER I, LNW, LP, LPB, LPL, LPP, NN ! ! Local parameters: ! ! I = DO-loop index ! LNW = LNEW-1 (output value of LNEW) ! LP = LIST pointer of the last neighbor of NB ! LPB = Pointer to NB as a neighbor of N0 ! LPL = Pointer to the last neighbor of N0 ! LPP = Pointer to the neighbor of N0 that precedes NB ! NN = Local copy of N ! NN = N ! ! Test for error 1. ! IF (N0 .LT. 1 .OR. N0 .GT. NN .OR. NB .LT. 1 .OR. & & NB .GT. NN .OR. NN .LT. 3) THEN LPH = -1 RETURN ENDIF ! ! Find pointers to neighbors of N0: ! ! LPL points to the last neighbor, ! LPP points to the neighbor NP preceding NB, and ! LPB points to NB. ! LPL = LEND(N0) LPP = LPL LPB = LPTR(LPP) 1 IF (LIST(LPB) .EQ. NB) GO TO 2 LPP = LPB LPB = LPTR(LPP) IF (LPB .NE. LPL) GO TO 1 ! ! Test for error 2 (NB not found). ! IF (ABS(LIST(LPB)) .NE. NB) THEN LPH = -2 RETURN ENDIF ! ! NB is the last neighbor of N0. Make NP the new last ! neighbor and, if NB is a boundary node, then make N0 ! a boundary node. ! LEND(N0) = LPP LP = LEND(NB) IF (LIST(LP) .LT. 0) LIST(LPP) = -LIST(LPP) GO TO 3 ! ! NB is not the last neighbor of N0. If NB is a boundary ! node and N0 is not, then make N0 a boundary node with ! last neighbor NP. ! 2 LP = LEND(NB) IF (LIST(LP) .LT. 0 .AND. LIST(LPL) .GT. 0) THEN LEND(N0) = LPP LIST(LPP) = -LIST(LPP) ENDIF ! ! Update LPTR so that the neighbor following NB now fol- ! lows NP, and fill in the hole at location LPB. ! 3 LPTR(LPP) = LPTR(LPB) LNW = LNEW-1 LIST(LPB) = LIST(LNW) LPTR(LPB) = LPTR(LNW) DO 4 I = NN,1,-1 IF (LEND(I) .EQ. LNW) THEN LEND(I) = LPB GO TO 5 ENDIF 4 CONTINUE ! 5 DO 6 I = 1,LNW-1 IF (LPTR(I) .EQ. LNW) THEN LPTR(I) = LPB ENDIF 6 CONTINUE ! ! No errors encountered. ! LNEW = LNW LPH = LPB RETURN END SUBROUTINE DELNOD (K, N,X,Y,Z,LIST,LPTR,LEND,LNEW,LWK, IWK, IER) INTEGER K, N, LIST(*), LPTR(*), LEND(*), LNEW, LWK, IWK(2,*), IER REAL X(*), Y(*), Z(*) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 11/30/99 ! ! This subroutine deletes node K (along with all arcs ! incident on node K) from a triangulation of N nodes on the ! unit sphere, and inserts arcs as necessary to produce a ! triangulation of the remaining N-1 nodes. If a Delaunay ! triangulation is input, a Delaunay triangulation will ! result, and thus, DELNOD reverses the effect of a call to ! Subroutine ADDNOD. ! ! ! On input: ! ! K = Index (for X, Y, and Z) of the node to be ! deleted. 1 .LE. K .LE. N. ! ! K is not altered by this routine. ! ! N = Number of nodes in the triangulation on input. ! N .GE. 4. Note that N will be decremented ! following the deletion. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes in the triangula- ! tion. ! ! LIST,LPTR,LEND,LNEW = Data structure defining the ! triangulation. Refer to Sub- ! routine TRMESH. ! ! LWK = Number of columns reserved for IWK. LWK must ! be at least NNB-3, where NNB is the number of ! neighbors of node K, including an extra ! pseudo-node if K is a boundary node. ! ! IWK = Integer work array dimensioned 2 by LWK (or ! array of length .GE. 2*LWK). ! ! On output: ! ! N = Number of nodes in the triangulation on output. ! The input value is decremented unless 1 .LE. IER ! .LE. 4. ! ! X,Y,Z = Updated arrays containing nodal coordinates ! (with elements K+1,...,N+1 shifted up one ! position, thus overwriting element K) unless ! 1 .LE. IER .LE. 4. ! ! LIST,LPTR,LEND,LNEW = Updated triangulation data ! structure reflecting the dele- ! tion unless 1 .LE. IER .LE. 4. ! Note that the data structure ! may have been altered if IER > ! 3. ! ! LWK = Number of IWK columns required unless IER = 1 ! or IER = 3. ! ! IWK = Indexes of the endpoints of the new arcs added ! unless LWK = 0 or 1 .LE. IER .LE. 4. (Arcs ! are associated with columns, or pairs of ! adjacent elements if IWK is declared as a ! singly-subscripted array.) ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if K or N is outside its valid range ! or LWK < 0 on input. ! IER = 2 if more space is required in IWK. ! Refer to LWK. ! IER = 3 if the triangulation data structure is ! invalid on input. ! IER = 4 if K indexes an interior node with ! four or more neighbors, none of which ! can be swapped out due to collineari- ! ty, and K cannot therefore be deleted. ! IER = 5 if an error flag (other than IER = 1) ! was returned by OPTIM. An error ! message is written to the standard ! output unit in this case. ! IER = 6 if error flag 1 was returned by OPTIM. ! This is not necessarily an error, but ! the arcs may not be optimal. ! ! Note that the deletion may result in all remaining nodes ! being collinear. This situation is not flagged. ! ! Modules required by DELNOD: DELNB, LEFT, LSTPTR, NBCNT, ! OPTIM, SWAP, SWPTST ! ! Intrinsic function called by DELNOD: ABS ! !*********************************************************** ! INTEGER LSTPTR, NBCNT INTEGER I, IERR, IWL, J, LNW, LP, LP21, LPF, LPH, LPL, & & LPL2, LPN, LWKL, N1, N2, NFRST, NIT, NL, NN, NNB, NR LOGICAL LEFT LOGICAL BDRY REAL X1, X2, XL, XR, Y1, Y2, YL, YR, Z1, Z2, ZL, ZR ! ! Local parameters: ! ! BDRY = Logical variable with value TRUE iff N1 is a ! boundary node ! I,J = DO-loop indexes ! IERR = Error flag returned by OPTIM ! IWL = Number of IWK columns containing arcs ! LNW = Local copy of LNEW ! LP = LIST pointer ! LP21 = LIST pointer returned by SWAP ! LPF,LPL = Pointers to the first and last neighbors of N1 ! LPH = Pointer (or flag) returned by DELNB ! LPL2 = Pointer to the last neighbor of N2 ! LPN = Pointer to a neighbor of N1 ! LWKL = Input value of LWK ! N1 = Local copy of K ! N2 = Neighbor of N1 ! NFRST = First neighbor of N1: LIST(LPF) ! NIT = Number of iterations in OPTIM ! NR,NL = Neighbors of N1 preceding (to the right of) and ! following (to the left of) N2, respectively ! NN = Number of nodes in the triangulation ! NNB = Number of neighbors of N1 (including a pseudo- ! node representing the boundary if N1 is a ! boundary node) ! X1,Y1,Z1 = Coordinates of N1 ! X2,Y2,Z2 = Coordinates of N2 ! XL,YL,ZL = Coordinates of NL ! XR,YR,ZR = Coordinates of NR ! ! ! Set N1 to K and NNB to the number of neighbors of N1 (plus ! one if N1 is a boundary node), and test for errors. LPF ! and LPL are LIST indexes of the first and last neighbors ! of N1, IWL is the number of IWK columns containing arcs, ! and BDRY is TRUE iff N1 is a boundary node. ! N1 = K NN = N IF (N1 .LT. 1 .OR. N1 .GT. NN .OR. NN .LT. 4 .OR. & & LWK .LT. 0) GO TO 21 LPL = LEND(N1) LPF = LPTR(LPL) NNB = NBCNT(LPL,LPTR) BDRY = LIST(LPL) .LT. 0 IF (BDRY) NNB = NNB + 1 IF (NNB .LT. 3) GO TO 23 LWKL = LWK LWK = NNB - 3 IF (LWKL .LT. LWK) GO TO 22 IWL = 0 IF (NNB .EQ. 3) GO TO 3 ! ! Initialize for loop on arcs N1-N2 for neighbors N2 of N1, ! beginning with the second neighbor. NR and NL are the ! neighbors preceding and following N2, respectively, and ! LP indexes NL. The loop is exited when all possible ! swaps have been applied to arcs incident on N1. ! X1 = X(N1) Y1 = Y(N1) Z1 = Z(N1) NFRST = LIST(LPF) NR = NFRST XR = X(NR) YR = Y(NR) ZR = Z(NR) LP = LPTR(LPF) N2 = LIST(LP) X2 = X(N2) Y2 = Y(N2) Z2 = Z(N2) LP = LPTR(LP) ! ! Top of loop: set NL to the neighbor following N2. ! 1 NL = ABS(LIST(LP)) IF (NL .EQ. NFRST .AND. BDRY) GO TO 3 XL = X(NL) YL = Y(NL) ZL = Z(NL) ! ! Test for a convex quadrilateral. To avoid an incorrect ! test caused by collinearity, use the fact that if N1 ! is a boundary node, then N1 LEFT NR->NL and if N2 is ! a boundary node, then N2 LEFT NL->NR. ! LPL2 = LEND(N2) IF ( .NOT. ((BDRY .OR. LEFT(XR,YR,ZR,XL,YL,ZL,X1,Y1, & & Z1)) .AND. (LIST(LPL2) .LT. 0 .OR. & & LEFT(XL,YL,ZL,XR,YR,ZR,X2,Y2,Z2))) ) THEN ! ! Nonconvex quadrilateral -- no swap is possible. ! NR = N2 XR = X2 YR = Y2 ZR = Z2 GO TO 2 ENDIF ! ! The quadrilateral defined by adjacent triangles ! (N1,N2,NL) and (N2,N1,NR) is convex. Swap in ! NL-NR and store it in IWK unless NL and NR are ! already adjacent, in which case the swap is not ! possible. Indexes larger than N1 must be decremented ! since N1 will be deleted from X, Y, and Z. ! CALL SWAP (NL,NR,N1,N2, LIST,LPTR,LEND, LP21) IF (LP21 .EQ. 0) THEN NR = N2 XR = X2 YR = Y2 ZR = Z2 GO TO 2 ENDIF IWL = IWL + 1 IF (NL .LE. N1) THEN IWK(1,IWL) = NL ELSE IWK(1,IWL) = NL - 1 ENDIF IF (NR .LE. N1) THEN IWK(2,IWL) = NR ELSE IWK(2,IWL) = NR - 1 ENDIF ! ! Recompute the LIST indexes and NFRST, and decrement NNB. ! LPL = LEND(N1) NNB = NNB - 1 IF (NNB .EQ. 3) GO TO 3 LPF = LPTR(LPL) NFRST = LIST(LPF) LP = LSTPTR(LPL,NL,LIST,LPTR) IF (NR .EQ. NFRST) GO TO 2 ! ! NR is not the first neighbor of N1. ! Back up and test N1-NR for a swap again: Set N2 to ! NR and NR to the previous neighbor of N1 -- the ! neighbor of NR which follows N1. LP21 points to NL ! as a neighbor of NR. ! N2 = NR X2 = XR Y2 = YR Z2 = ZR LP21 = LPTR(LP21) LP21 = LPTR(LP21) NR = ABS(LIST(LP21)) XR = X(NR) YR = Y(NR) ZR = Z(NR) GO TO 1 ! ! Bottom of loop -- test for termination of loop. ! 2 IF (N2 .EQ. NFRST) GO TO 3 N2 = NL X2 = XL Y2 = YL Z2 = ZL LP = LPTR(LP) GO TO 1 ! ! Delete N1 and all its incident arcs. If N1 is an interior ! node and either NNB > 3 or NNB = 3 and N2 LEFT NR->NL, ! then N1 must be separated from its neighbors by a plane ! containing the origin -- its removal reverses the effect ! of a call to COVSPH, and all its neighbors become ! boundary nodes. This is achieved by treating it as if ! it were a boundary node (setting BDRY to TRUE, changing ! a sign in LIST, and incrementing NNB). ! 3 IF (.NOT. BDRY) THEN IF (NNB .GT. 3) THEN BDRY = .TRUE. ELSE LPF = LPTR(LPL) NR = LIST(LPF) LP = LPTR(LPF) N2 = LIST(LP) NL = LIST(LPL) BDRY = LEFT(X(NR),Y(NR),Z(NR),X(NL),Y(NL),Z(NL), & & X(N2),Y(N2),Z(N2)) ENDIF IF (BDRY) THEN ! ! IF a boundary node already exists, then N1 and its ! neighbors cannot be converted to boundary nodes. ! (They must be collinear.) This is a problem if ! NNB > 3. ! DO 4 I = 1,NN IF (LIST(LEND(I)) .LT. 0) THEN BDRY = .FALSE. GO TO 5 ENDIF 4 CONTINUE LIST(LPL) = -LIST(LPL) NNB = NNB + 1 ENDIF ENDIF 5 IF (.NOT. BDRY .AND. NNB .GT. 3) GO TO 24 ! ! Initialize for loop on neighbors. LPL points to the last ! neighbor of N1. LNEW is stored in local variable LNW. ! LP = LPL LNW = LNEW ! ! Loop on neighbors N2 of N1, beginning with the first. ! 6 LP = LPTR(LP) N2 = ABS(LIST(LP)) CALL DELNB (N2,N1,N, LIST,LPTR,LEND,LNW, LPH) IF (LPH .LT. 0) GO TO 23 ! ! LP and LPL may require alteration. ! IF (LPL .EQ. LNW) LPL = LPH IF (LP .EQ. LNW) LP = LPH IF (LP .NE. LPL) GO TO 6 ! ! Delete N1 from X, Y, Z, and LEND, and remove its adjacency ! list from LIST and LPTR. LIST entries (nodal indexes) ! which are larger than N1 must be decremented. ! NN = NN - 1 IF (N1 .GT. NN) GO TO 9 DO 7 I = N1,NN X(I) = X(I+1) Y(I) = Y(I+1) Z(I) = Z(I+1) LEND(I) = LEND(I+1) 7 CONTINUE ! DO 8 I = 1,LNW-1 IF (LIST(I) .GT. N1) LIST(I) = LIST(I) - 1 IF (LIST(I) .LT. -N1) LIST(I) = LIST(I) + 1 8 CONTINUE ! ! For LPN = first to last neighbors of N1, delete the ! preceding neighbor (indexed by LP). ! ! Each empty LIST,LPTR location LP is filled in with the ! values at LNW-1, and LNW is decremented. All pointers ! (including those in LPTR and LEND) with value LNW-1 ! must be changed to LP. ! ! LPL points to the last neighbor of N1. ! 9 IF (BDRY) NNB = NNB - 1 LPN = LPL DO 13 J = 1,NNB LNW = LNW - 1 LP = LPN LPN = LPTR(LP) LIST(LP) = LIST(LNW) LPTR(LP) = LPTR(LNW) IF (LPTR(LPN) .EQ. LNW) LPTR(LPN) = LP IF (LPN .EQ. LNW) LPN = LP DO 10 I = NN,1,-1 IF (LEND(I) .EQ. LNW) THEN LEND(I) = LP GO TO 11 ENDIF 10 CONTINUE ! 11 DO 12 I = LNW-1,1,-1 IF (LPTR(I) .EQ. LNW) LPTR(I) = LP 12 CONTINUE 13 CONTINUE ! ! Update N and LNEW, and optimize the patch of triangles ! containing K (on input) by applying swaps to the arcs ! in IWK. ! N = NN LNEW = LNW IF (IWL .GT. 0) THEN NIT = 4*IWL CALL OPTIM (X,Y,Z,IWL, LIST,LPTR,LEND,NIT,IWK, IERR) IF (IERR .NE. 0 .AND. IERR .NE. 1) GO TO 25 IF (IERR .EQ. 1) GO TO 26 ENDIF ! ! Successful termination. ! IER = 0 RETURN ! ! Invalid input parameter. ! 21 IER = 1 RETURN ! ! Insufficient space reserved for IWK. ! 22 IER = 2 RETURN ! ! Invalid triangulation data structure. NNB < 3 on input or ! N2 is a neighbor of N1 but N1 is not a neighbor of N2. ! 23 IER = 3 RETURN ! ! N1 is interior but NNB could not be reduced to 3. ! 24 IER = 4 RETURN ! ! Error flag (other than 1) returned by OPTIM. ! 25 IER = 5 WRITE (*,100) NIT, IERR 100 FORMAT (//5X,'*** Error in OPTIM (called from ', & & 'DELNOD): NIT = ',I4,', IER = ',I1,' ***'/) RETURN ! ! Error flag 1 returned by OPTIM. ! 26 IER = 6 RETURN END SUBROUTINE EDGE (IN1,IN2,X,Y,Z, LWK,IWK,LIST,LPTR, LEND, IER) INTEGER IN1, IN2, LWK, IWK(2,*), LIST(*), LPTR(*), LEND(*), IER REAL X(*), Y(*), Z(*) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/30/98 ! ! Given a triangulation of N nodes and a pair of nodal ! indexes IN1 and IN2, this routine swaps arcs as necessary ! to force IN1 and IN2 to be adjacent. Only arcs which ! intersect IN1-IN2 are swapped out. If a Delaunay triangu- ! lation is input, the resulting triangulation is as close ! as possible to a Delaunay triangulation in the sense that ! all arcs other than IN1-IN2 are locally optimal. ! ! A sequence of calls to EDGE may be used to force the ! presence of a set of edges defining the boundary of a non- ! convex and/or multiply connected region, or to introduce ! barriers into the triangulation. Note that Subroutine ! GETNP will not necessarily return closest nodes if the ! triangulation has been constrained by a call to EDGE. ! However, this is appropriate in some applications, such ! as triangle-based interpolation on a nonconvex domain. ! ! ! On input: ! ! IN1,IN2 = Indexes (of X, Y, and Z) in the range 1 to ! N defining a pair of nodes to be connected ! by an arc. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes. ! ! The above parameters are not altered by this routine. ! ! LWK = Number of columns reserved for IWK. This must ! be at least NI -- the number of arcs that ! intersect IN1-IN2. (NI is bounded by N-3.) ! ! IWK = Integer work array of length at least 2*LWK. ! ! LIST,LPTR,LEND = Data structure defining the trian- ! gulation. Refer to Subroutine ! TRMESH. ! ! On output: ! ! LWK = Number of arcs which intersect IN1-IN2 (but ! not more than the input value of LWK) unless ! IER = 1 or IER = 3. LWK = 0 if and only if ! IN1 and IN2 were adjacent (or LWK=0) on input. ! ! IWK = Array containing the indexes of the endpoints ! of the new arcs other than IN1-IN2 unless ! IER > 0 or LWK = 0. New arcs to the left of ! IN1->IN2 are stored in the first K-1 columns ! (left portion of IWK), column K contains ! zeros, and new arcs to the right of IN1->IN2 ! occupy columns K+1,...,LWK. (K can be deter- ! mined by searching IWK for the zeros.) ! ! LIST,LPTR,LEND = Data structure updated if necessary ! to reflect the presence of an arc ! connecting IN1 and IN2 unless IER > ! 0. The data structure has been ! altered if IER >= 4. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if IN1 < 1, IN2 < 1, IN1 = IN2, ! or LWK < 0 on input. ! IER = 2 if more space is required in IWK. ! Refer to LWK. ! IER = 3 if IN1 and IN2 could not be connected ! due to either an invalid data struc- ! ture or collinear nodes (and floating ! point error). ! IER = 4 if an error flag other than IER = 1 ! was returned by OPTIM. ! IER = 5 if error flag 1 was returned by OPTIM. ! This is not necessarily an error, but ! the arcs other than IN1-IN2 may not ! be optimal. ! ! An error message is written to the standard output unit ! in the case of IER = 3 or IER = 4. ! ! Modules required by EDGE: LEFT, LSTPTR, OPTIM, SWAP, ! SWPTST ! ! Intrinsic function called by EDGE: ABS ! !*********************************************************** ! LOGICAL LEFT INTEGER I, IERR, IWC, IWCP1, IWEND, IWF, IWL, LFT, LP, & & LP21, LPL, N0, N1, N1FRST, N1LST, N2, NEXT, & & NIT, NL, NR REAL DP12, DP1L, DP1R, DP2L, DP2R, X0, X1, X2, Y0, & & Y1, Y2, Z0, Z1, Z2 ! ! Local parameters: ! ! DPij = Dot product ! I = DO-loop index and column index for IWK ! IERR = Error flag returned by Subroutine OPTIM ! IWC = IWK index between IWF and IWL -- NL->NR is ! stored in IWK(1,IWC)->IWK(2,IWC) ! IWCP1 = IWC + 1 ! IWEND = Input or output value of LWK ! IWF = IWK (column) index of the first (leftmost) arc ! which intersects IN1->IN2 ! IWL = IWK (column) index of the last (rightmost) are ! which intersects IN1->IN2 ! LFT = Flag used to determine if a swap results in the ! new arc intersecting IN1-IN2 -- LFT = 0 iff ! N0 = IN1, LFT = -1 implies N0 LEFT IN1->IN2, ! and LFT = 1 implies N0 LEFT IN2->IN1 ! LP = List pointer (index for LIST and LPTR) ! LP21 = Unused parameter returned by SWAP ! LPL = Pointer to the last neighbor of IN1 or NL ! N0 = Neighbor of N1 or node opposite NR->NL ! N1,N2 = Local copies of IN1 and IN2 ! N1FRST = First neighbor of IN1 ! N1LST = (Signed) last neighbor of IN1 ! NEXT = Node opposite NL->NR ! NIT = Flag or number of iterations employed by OPTIM ! NL,NR = Endpoints of an arc which intersects IN1-IN2 ! with NL LEFT IN1->IN2 ! X0,Y0,Z0 = Coordinates of N0 ! X1,Y1,Z1 = Coordinates of IN1 ! X2,Y2,Z2 = Coordinates of IN2 ! ! ! Store IN1, IN2, and LWK in local variables and test for ! errors. ! N1 = IN1 N2 = IN2 IWEND = LWK IF (N1 .LT. 1 .OR. N2 .LT. 1 .OR. N1 .EQ. N2 .OR. & & IWEND .LT. 0) GO TO 31 ! ! Test for N2 as a neighbor of N1. LPL points to the last ! neighbor of N1. ! LPL = LEND(N1) N0 = ABS(LIST(LPL)) LP = LPL 1 IF (N0 .EQ. N2) GO TO 30 LP = LPTR(LP) N0 = LIST(LP) IF (LP .NE. LPL) GO TO 1 ! ! Initialize parameters. ! IWL = 0 NIT = 0 ! ! Store the coordinates of N1 and N2. ! 2 X1 = X(N1) Y1 = Y(N1) Z1 = Z(N1) X2 = X(N2) Y2 = Y(N2) Z2 = Z(N2) ! ! Set NR and NL to adjacent neighbors of N1 such that ! NR LEFT N2->N1 and NL LEFT N1->N2, ! (NR Forward N1->N2 or NL Forward N1->N2), and ! (NR Forward N2->N1 or NL Forward N2->N1). ! ! Initialization: Set N1FRST and N1LST to the first and ! (signed) last neighbors of N1, respectively, and ! initialize NL to N1FRST. ! LPL = LEND(N1) N1LST = LIST(LPL) LP = LPTR(LPL) N1FRST = LIST(LP) NL = N1FRST IF (N1LST .LT. 0) GO TO 4 ! ! N1 is an interior node. Set NL to the first candidate ! for NR (NL LEFT N2->N1). ! 3 IF (LEFT(X2,Y2,Z2,X1,Y1,Z1,X(NL),Y(NL),Z(NL))) GO TO 4 LP = LPTR(LP) NL = LIST(LP) IF (NL .NE. N1FRST) GO TO 3 ! ! All neighbors of N1 are strictly left of N1->N2. ! GO TO 5 ! ! NL = LIST(LP) LEFT N2->N1. Set NR to NL and NL to the ! following neighbor of N1. ! 4 NR = NL LP = LPTR(LP) NL = ABS(LIST(LP)) IF (LEFT(X1,Y1,Z1,X2,Y2,Z2,X(NL),Y(NL),Z(NL)) ) THEN ! ! NL LEFT N1->N2 and NR LEFT N2->N1. The Forward tests ! are employed to avoid an error associated with ! collinear nodes. ! DP12 = X1*X2 + Y1*Y2 + Z1*Z2 DP1L = X1*X(NL) + Y1*Y(NL) + Z1*Z(NL) DP2L = X2*X(NL) + Y2*Y(NL) + Z2*Z(NL) DP1R = X1*X(NR) + Y1*Y(NR) + Z1*Z(NR) DP2R = X2*X(NR) + Y2*Y(NR) + Z2*Z(NR) IF ( (DP2L-DP12*DP1L .GE. 0. .OR. & & DP2R-DP12*DP1R .GE. 0.) .AND. & & (DP1L-DP12*DP2L .GE. 0. .OR. & & DP1R-DP12*DP2R .GE. 0.) ) GO TO 6 ! ! NL-NR does not intersect N1-N2. However, there is ! another candidate for the first arc if NL lies on ! the line N1-N2. ! IF ( .NOT. LEFT(X2,Y2,Z2,X1,Y1,Z1,X(NL),Y(NL), & & Z(NL)) ) GO TO 5 ENDIF ! ! Bottom of loop. ! IF (NL .NE. N1FRST) GO TO 4 ! ! Either the triangulation is invalid or N1-N2 lies on the ! convex hull boundary and an edge NR->NL (opposite N1 and ! intersecting N1-N2) was not found due to floating point ! error. Try interchanging N1 and N2 -- NIT > 0 iff this ! has already been done. ! 5 IF (NIT .GT. 0) GO TO 33 NIT = 1 N1 = N2 N2 = IN1 GO TO 2 ! ! Store the ordered sequence of intersecting edges NL->NR in ! IWK(1,IWL)->IWK(2,IWL). ! 6 IWL = IWL + 1 IF (IWL .GT. IWEND) GO TO 32 IWK(1,IWL) = NL IWK(2,IWL) = NR ! ! Set NEXT to the neighbor of NL which follows NR. ! LPL = LEND(NL) LP = LPTR(LPL) ! ! Find NR as a neighbor of NL. The search begins with ! the first neighbor. ! 7 IF (LIST(LP) .EQ. NR) GO TO 8 LP = LPTR(LP) IF (LP .NE. LPL) GO TO 7 ! ! NR must be the last neighbor, and NL->NR cannot be a ! boundary edge. ! IF (LIST(LP) .NE. NR) GO TO 33 ! ! Set NEXT to the neighbor following NR, and test for ! termination of the store loop. ! 8 LP = LPTR(LP) NEXT = ABS(LIST(LP)) IF (NEXT .EQ. N2) GO TO 9 ! ! Set NL or NR to NEXT. ! IF ( LEFT(X1,Y1,Z1,X2,Y2,Z2,X(NEXT),Y(NEXT),Z(NEXT)) ) THEN NL = NEXT ELSE NR = NEXT ENDIF GO TO 6 ! ! IWL is the number of arcs which intersect N1-N2. ! Store LWK. ! 9 LWK = IWL IWEND = IWL ! ! Initialize for edge swapping loop -- all possible swaps ! are applied (even if the new arc again intersects ! N1-N2), arcs to the left of N1->N2 are stored in the ! left portion of IWK, and arcs to the right are stored in ! the right portion. IWF and IWL index the first and last ! intersecting arcs. ! IWF = 1 ! ! Top of loop -- set N0 to N1 and NL->NR to the first edge. ! IWC points to the arc currently being processed. LFT ! .LE. 0 iff N0 LEFT N1->N2. ! 10 LFT = 0 N0 = N1 X0 = X1 Y0 = Y1 Z0 = Z1 NL = IWK(1,IWF) NR = IWK(2,IWF) IWC = IWF ! ! Set NEXT to the node opposite NL->NR unless IWC is the ! last arc. ! 11 IF (IWC .EQ. IWL) GO TO 21 IWCP1 = IWC + 1 NEXT = IWK(1,IWCP1) IF (NEXT .NE. NL) GO TO 16 NEXT = IWK(2,IWCP1) ! ! NEXT RIGHT N1->N2 and IWC .LT. IWL. Test for a possible ! swap. ! IF ( .NOT. LEFT(X0,Y0,Z0,X(NR),Y(NR),Z(NR),X(NEXT), & & Y(NEXT),Z(NEXT)) ) GO TO 14 IF (LFT .GE. 0) GO TO 12 IF ( .NOT. LEFT(X(NL),Y(NL),Z(NL),X0,Y0,Z0,X(NEXT), & & Y(NEXT),Z(NEXT)) ) GO TO 14 ! ! Replace NL->NR with N0->NEXT. ! CALL SWAP (NEXT,N0,NL,NR, LIST,LPTR,LEND, LP21) IWK(1,IWC) = N0 IWK(2,IWC) = NEXT GO TO 15 ! ! Swap NL-NR for N0-NEXT, shift columns IWC+1,...,IWL to ! the left, and store N0-NEXT in the right portion of ! IWK. ! 12 CALL SWAP (NEXT,N0,NL,NR, LIST,LPTR,LEND, LP21) DO 13 I = IWCP1,IWL IWK(1,I-1) = IWK(1,I) IWK(2,I-1) = IWK(2,I) 13 CONTINUE IWK(1,IWL) = N0 IWK(2,IWL) = NEXT IWL = IWL - 1 NR = NEXT GO TO 11 ! ! A swap is not possible. Set N0 to NR. ! 14 N0 = NR X0 = X(N0) Y0 = Y(N0) Z0 = Z(N0) LFT = 1 ! ! Advance to the next arc. ! 15 NR = NEXT IWC = IWC + 1 GO TO 11 ! ! NEXT LEFT N1->N2, NEXT .NE. N2, and IWC .LT. IWL. ! Test for a possible swap. ! 16 IF ( .NOT. LEFT(X(NL),Y(NL),Z(NL),X0,Y0,Z0,X(NEXT), & & Y(NEXT),Z(NEXT)) ) GO TO 19 IF (LFT .LE. 0) GO TO 17 IF ( .NOT. LEFT(X0,Y0,Z0,X(NR),Y(NR),Z(NR),X(NEXT), & & Y(NEXT),Z(NEXT)) ) GO TO 19 ! ! Replace NL->NR with NEXT->N0. ! CALL SWAP (NEXT,N0,NL,NR, LIST,LPTR,LEND, LP21) IWK(1,IWC) = NEXT IWK(2,IWC) = N0 GO TO 20 ! ! Swap NL-NR for N0-NEXT, shift columns IWF,...,IWC-1 to ! the right, and store N0-NEXT in the left portion of ! IWK. ! 17 CALL SWAP (NEXT,N0,NL,NR, LIST,LPTR,LEND, LP21) DO 18 I = IWC-1,IWF,-1 IWK(1,I+1) = IWK(1,I) IWK(2,I+1) = IWK(2,I) 18 CONTINUE IWK(1,IWF) = N0 IWK(2,IWF) = NEXT IWF = IWF + 1 GO TO 20 ! ! A swap is not possible. Set N0 to NL. ! 19 N0 = NL X0 = X(N0) Y0 = Y(N0) Z0 = Z(N0) LFT = -1 ! ! Advance to the next arc. ! 20 NL = NEXT IWC = IWC + 1 GO TO 11 ! ! N2 is opposite NL->NR (IWC = IWL). ! 21 IF (N0 .EQ. N1) GO TO 24 IF (LFT .LT. 0) GO TO 22 ! ! N0 RIGHT N1->N2. Test for a possible swap. ! IF ( .NOT. LEFT(X0,Y0,Z0,X(NR),Y(NR),Z(NR),X2,Y2,Z2) ) & & GO TO 10 ! ! Swap NL-NR for N0-N2 and store N0-N2 in the right ! portion of IWK. ! CALL SWAP (N2,N0,NL,NR, LIST,LPTR,LEND, LP21) IWK(1,IWL) = N0 IWK(2,IWL) = N2 IWL = IWL - 1 GO TO 10 ! ! N0 LEFT N1->N2. Test for a possible swap. ! 22 IF ( .NOT. LEFT(X(NL),Y(NL),Z(NL),X0,Y0,Z0,X2,Y2,Z2) ) & & GO TO 10 ! ! Swap NL-NR for N0-N2, shift columns IWF,...,IWL-1 to the ! right, and store N0-N2 in the left portion of IWK. ! CALL SWAP (N2,N0,NL,NR, LIST,LPTR,LEND, LP21) I = IWL 23 IWK(1,I) = IWK(1,I-1) IWK(2,I) = IWK(2,I-1) I = I - 1 IF (I .GT. IWF) GO TO 23 IWK(1,IWF) = N0 IWK(2,IWF) = N2 IWF = IWF + 1 GO TO 10 ! ! IWF = IWC = IWL. Swap out the last arc for N1-N2 and ! store zeros in IWK. ! 24 CALL SWAP (N2,N1,NL,NR, LIST,LPTR,LEND, LP21) IWK(1,IWC) = 0 IWK(2,IWC) = 0 ! ! Optimization procedure -- ! IER = 0 IF (IWC .GT. 1) THEN ! ! Optimize the set of new arcs to the left of IN1->IN2. ! NIT = 4*(IWC-1) CALL OPTIM (X,Y,Z,IWC-1, LIST,LPTR,LEND,NIT, IWK, IERR) IF (IERR .NE. 0 .AND. IERR .NE. 1) GO TO 34 IF (IERR .EQ. 1) IER = 5 ENDIF IF (IWC .LT. IWEND) THEN ! ! Optimize the set of new arcs to the right of IN1->IN2. ! NIT = 4*(IWEND-IWC) CALL OPTIM (X,Y,Z,IWEND-IWC, LIST,LPTR,LEND,NIT, & & IWK(1,IWC+1), IERR) IF (IERR .NE. 0 .AND. IERR .NE. 1) GO TO 34 IF (IERR .EQ. 1) GO TO 35 ENDIF IF (IER .EQ. 5) GO TO 35 ! ! Successful termination (IER = 0). ! RETURN ! ! IN1 and IN2 were adjacent on input. ! 30 IER = 0 RETURN ! ! Invalid input parameter. ! 31 IER = 1 RETURN ! ! Insufficient space reserved for IWK. ! 32 IER = 2 RETURN ! ! Invalid triangulation data structure or collinear nodes ! on convex hull boundary. ! 33 IER = 3 WRITE (*,130) IN1, IN2 130 FORMAT (//5X,'*** Error in EDGE: Invalid triangula', & & 'tion or null triangles on boundary'/ & & 9X,'IN1 =',I4,', IN2=',I4/) RETURN ! ! Error flag (other than 1) returned by OPTIM. ! 34 IER = 4 WRITE (*,140) NIT, IERR 140 FORMAT (//5X,'*** Error in OPTIM (called from EDGE):', & & ' NIT = ',I4,', IER = ',I1,' ***'/) RETURN ! ! Error flag 1 returned by OPTIM. ! 35 IER = 5 RETURN END SUBROUTINE GETNP (X,Y,Z,LIST,LPTR,LEND,L, NPTS, DF, IER) INTEGER LIST(*), LPTR(*), LEND(*), L, NPTS(L), IER REAL X(*), Y(*), Z(*), DF ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/28/98 ! ! Given a Delaunay triangulation of N nodes on the unit ! sphere and an array NPTS containing the indexes of L-1 ! nodes ordered by angular distance from NPTS(1), this sub- ! routine sets NPTS(L) to the index of the next node in the ! sequence -- the node, other than NPTS(1),...,NPTS(L-1), ! that is closest to NPTS(1). Thus, the ordered sequence ! of K closest nodes to N1 (including N1) may be determined ! by K-1 calls to GETNP with NPTS(1) = N1 and L = 2,3,...,K ! for K .GE. 2. ! ! The algorithm uses the property of a Delaunay triangula- ! tion that the K-th closest node to N1 is a neighbor of one ! of the K-1 closest nodes to N1. ! ! ! On input: ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes. ! ! LIST,LPTR,LEND = Triangulation data structure. Re- ! fer to Subroutine TRMESH. ! ! L = Number of nodes in the sequence on output. 2 ! .LE. L .LE. N. ! ! The above parameters are not altered by this routine. ! ! NPTS = Array of length .GE. L containing the indexes ! of the L-1 closest nodes to NPTS(1) in the ! first L-1 locations. ! ! On output: ! ! NPTS = Array updated with the index of the L-th ! closest node to NPTS(1) in position L unless ! IER = 1. ! ! DF = Value of an increasing function (negative cos- ! ine) of the angular distance between NPTS(1) ! and NPTS(L) unless IER = 1. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if L < 2. ! ! Modules required by GETNP: None ! ! Intrinsic function called by GETNP: ABS ! !*********************************************************** ! INTEGER I, LM1, LP, LPL, N1, NB, NI, NP REAL DNB, DNP, X1, Y1, Z1 ! ! Local parameters: ! ! DNB,DNP = Negative cosines of the angular distances from ! N1 to NB and to NP, respectively ! I = NPTS index and DO-loop index ! LM1 = L-1 ! LP = LIST pointer of a neighbor of NI ! LPL = Pointer to the last neighbor of NI ! N1 = NPTS(1) ! NB = Neighbor of NI and candidate for NP ! NI = NPTS(I) ! NP = Candidate for NPTS(L) ! X1,Y1,Z1 = Coordinates of N1 ! LM1 = L - 1 IF (LM1 .LT. 1) GO TO 6 IER = 0 ! ! Store N1 = NPTS(1) and mark the elements of NPTS. ! N1 = NPTS(1) X1 = X(N1) Y1 = Y(N1) Z1 = Z(N1) DO 1 I = 1,LM1 NI = NPTS(I) LEND(NI) = -LEND(NI) 1 CONTINUE ! ! Candidates for NP = NPTS(L) are the unmarked neighbors ! of nodes in NPTS. DNP is initially greater than -cos(PI) ! (the maximum distance). ! DNP = 2. ! ! Loop on nodes NI in NPTS. ! DO 4 I = 1,LM1 NI = NPTS(I) LPL = -LEND(NI) LP = LPL ! ! Loop on neighbors NB of NI. ! 2 NB = ABS(LIST(LP)) IF (LEND(NB) .LT. 0) GO TO 3 ! ! NB is an unmarked neighbor of NI. Replace NP if NB is ! closer to N1. ! DNB = -(X(NB)*X1 + Y(NB)*Y1 + Z(NB)*Z1) IF (DNB .GE. DNP) GO TO 3 NP = NB DNP = DNB 3 LP = LPTR(LP) IF (LP .NE. LPL) GO TO 2 4 CONTINUE NPTS(L) = NP DF = DNP ! ! Unmark the elements of NPTS. ! DO 5 I = 1,LM1 NI = NPTS(I) LEND(NI) = -LEND(NI) 5 CONTINUE RETURN ! ! L is outside its valid range. ! 6 IER = 1 RETURN END SUBROUTINE INSERT (K,LP, LIST,LPTR,LNEW ) INTEGER K, LP, LIST(*), LPTR(*), LNEW ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/17/96 ! ! This subroutine inserts K as a neighbor of N1 following ! N2, where LP is the LIST pointer of N2 as a neighbor of ! N1. Note that, if N2 is the last neighbor of N1, K will ! become the first neighbor (even if N1 is a boundary node). ! ! This routine is identical to the similarly named routine ! in TRIPACK. ! ! ! On input: ! ! K = Index of the node to be inserted. ! ! LP = LIST pointer of N2 as a neighbor of N1. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LNEW = Data structure defining the trian- ! gulation. Refer to Subroutine ! TRMESH. ! ! On output: ! ! LIST,LPTR,LNEW = Data structure updated with the ! addition of node K. ! ! Modules required by INSERT: None ! !*********************************************************** ! INTEGER LSAV ! LSAV = LPTR(LP) LPTR(LP) = LNEW LIST(LNEW) = K LPTR(LNEW) = LSAV LNEW = LNEW + 1 RETURN END LOGICAL FUNCTION INSIDE (P,LV,XV,YV,ZV,NV,LISTV, IER) INTEGER LV, NV, LISTV(NV), IER REAL P(3), XV(LV), YV(LV), ZV(LV) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 12/27/93 ! ! This function locates a point P relative to a polygonal ! region R on the surface of the unit sphere, returning ! INSIDE = TRUE if and only if P is contained in R. R is ! defined by a cyclically ordered sequence of vertices which ! form a positively-oriented simple closed curve. Adjacent ! vertices need not be distinct but the curve must not be ! self-intersecting. Also, while polygon edges are by defi- ! nition restricted to a single hemisphere, R is not so ! restricted. Its interior is the region to the left as the ! vertices are traversed in order. ! ! The algorithm consists of selecting a point Q in R and ! then finding all points at which the great circle defined ! by P and Q intersects the boundary of R. P lies inside R ! if and only if there is an even number of intersection ! points between Q and P. Q is taken to be a point immedi- ! ately to the left of a directed boundary edge -- the first ! one that results in no consistency-check failures. ! ! If P is close to the polygon boundary, the problem is ! ill-conditioned and the decision may be incorrect. Also, ! an incorrect decision may result from a poor choice of Q ! (if, for example, a boundary edge lies on the great cir- ! cle defined by P and Q). A more reliable result could be ! obtained by a sequence of calls to INSIDE with the ver- ! tices cyclically permuted before each call (to alter the ! choice of Q). ! ! ! On input: ! ! P = Array of length 3 containing the Cartesian ! coordinates of the point (unit vector) to be ! located. ! ! LV = Length of arrays XV, YV, and ZV. ! ! XV,YV,ZV = Arrays of length LV containing the Carte- ! sian coordinates of unit vectors (points ! on the unit sphere). These values are ! not tested for validity. ! ! NV = Number of vertices in the polygon. 3 .LE. NV ! .LE. LV. ! ! LISTV = Array of length NV containing the indexes ! (for XV, YV, and ZV) of a cyclically-ordered ! (and CCW-ordered) sequence of vertices that ! define R. The last vertex (indexed by ! LISTV(NV)) is followed by the first (indexed ! by LISTV(1)). LISTV entries must be in the ! range 1 to LV. ! ! Input parameters are not altered by this function. ! ! On output: ! ! INSIDE = TRUE if and only if P lies inside R unless ! IER .NE. 0, in which case the value is not ! altered. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if LV or NV is outside its valid ! range. ! IER = 2 if a LISTV entry is outside its valid ! range. ! IER = 3 if the polygon boundary was found to ! be self-intersecting. This error will ! not necessarily be detected. ! IER = 4 if every choice of Q (one for each ! boundary edge) led to failure of some ! internal consistency check. The most ! likely cause of this error is invalid ! input: P = (0,0,0), a null or self- ! intersecting polygon, etc. ! ! Module required by INSIDE: INTRSC ! ! Intrinsic function called by INSIDE: SQRT ! !*********************************************************** ! INTEGER I1, I2, IERR, IMX, K, K0, N, NI LOGICAL EVEN, LFT1, LFT2, PINR, QINR REAL B(3), BP, BQ, CN(3), D, EPS, PN(3), Q(3), & & QN(3), QNRM, V1(3), V2(3), VN(3), VNRM ! ! Local parameters: ! ! B = Intersection point between the boundary and ! the great circle defined by P and Q ! BP,BQ = and , respectively, maximized over ! intersection points B that lie between P and ! Q (on the shorter arc) -- used to find the ! closest intersection points to P and Q ! CN = Q X P = normal to the plane of P and Q ! D = Dot product or ! EPS = Parameter used to define Q as the point whose ! orthogonal distance to (the midpoint of) ! boundary edge V1->V2 is approximately EPS/ ! (2*Cos(A/2)), where = Cos(A). ! EVEN = TRUE iff an even number of intersection points ! lie between P and Q (on the shorter arc) ! I1,I2 = Indexes (LISTV elements) of a pair of adjacent ! boundary vertices (endpoints of a boundary ! edge) ! IERR = Error flag for calls to INTRSC (not tested) ! IMX = Local copy of LV and maximum value of I1 and ! I2 ! K = DO-loop index and LISTV index ! K0 = LISTV index of the first endpoint of the ! boundary edge used to compute Q ! LFT1,LFT2 = Logical variables associated with I1 and I2 in ! the boundary traversal: TRUE iff the vertex ! is strictly to the left of Q->P ( > 0) ! N = Local copy of NV ! NI = Number of intersections (between the boundary ! curve and the great circle P-Q) encountered ! PINR = TRUE iff P is to the left of the directed ! boundary edge associated with the closest ! intersection point to P that lies between P ! and Q (a left-to-right intersection as ! viewed from Q), or there is no intersection ! between P and Q (on the shorter arc) ! PN,QN = P X CN and CN X Q, respectively: used to ! locate intersections B relative to arc Q->P ! Q = (V1 + V2 + EPS*VN/VNRM)/QNRM, where V1->V2 is ! the boundary edge indexed by LISTV(K0) -> ! LISTV(K0+1) ! QINR = TRUE iff Q is to the left of the directed ! boundary edge associated with the closest ! intersection point to Q that lies between P ! and Q (a right-to-left intersection as ! viewed from Q), or there is no intersection ! between P and Q (on the shorter arc) ! QNRM = Euclidean norm of V1+V2+EPS*VN/VNRM used to ! compute (normalize) Q ! V1,V2 = Vertices indexed by I1 and I2 in the boundary ! traversal ! VN = V1 X V2, where V1->V2 is the boundary edge ! indexed by LISTV(K0) -> LISTV(K0+1) ! VNRM = Euclidean norm of VN ! ! DATA EPS/1.E-4/ DATA EPS/0.5E-4/ ! ! Store local parameters, test for error 1, and initialize ! K0. ! IMX = LV N = NV IF (N .LT. 3 .OR. N .GT. IMX) GO TO 11 K0 = 0 I1 = LISTV(1) IF (I1 .LT. 1 .OR. I1 .GT. IMX) GO TO 12 ! ! Increment K0 and set Q to a point immediately to the left ! of the midpoint of edge V1->V2 = LISTV(K0)->LISTV(K0+1): ! Q = (V1 + V2 + EPS*VN/VNRM)/QNRM, where VN = V1 X V2. ! 1 K0 = K0 + 1 IF (K0 .GT. N) GO TO 14 I1 = LISTV(K0) IF (K0 .LT. N) THEN I2 = LISTV(K0+1) ELSE I2 = LISTV(1) ENDIF IF (I2 .LT. 1 .OR. I2 .GT. IMX) GO TO 12 VN(1) = YV(I1)*ZV(I2) - ZV(I1)*YV(I2) VN(2) = ZV(I1)*XV(I2) - XV(I1)*ZV(I2) VN(3) = XV(I1)*YV(I2) - YV(I1)*XV(I2) VNRM = SQRT(VN(1)*VN(1) + VN(2)*VN(2) + VN(3)*VN(3)) IF (VNRM .EQ. 0.) GO TO 1 Q(1) = XV(I1) + XV(I2) + EPS*VN(1)/VNRM Q(2) = YV(I1) + YV(I2) + EPS*VN(2)/VNRM Q(3) = ZV(I1) + ZV(I2) + EPS*VN(3)/VNRM QNRM = SQRT(Q(1)*Q(1) + Q(2)*Q(2) + Q(3)*Q(3)) Q(1) = Q(1)/QNRM Q(2) = Q(2)/QNRM Q(3) = Q(3)/QNRM ! ! Compute CN = Q X P, PN = P X CN, and QN = CN X Q. ! CN(1) = Q(2)*P(3) - Q(3)*P(2) CN(2) = Q(3)*P(1) - Q(1)*P(3) CN(3) = Q(1)*P(2) - Q(2)*P(1) IF (CN(1) .EQ. 0. .AND. CN(2) .EQ. 0. .AND. & & CN(3) .EQ. 0.) GO TO 1 PN(1) = P(2)*CN(3) - P(3)*CN(2) PN(2) = P(3)*CN(1) - P(1)*CN(3) PN(3) = P(1)*CN(2) - P(2)*CN(1) QN(1) = CN(2)*Q(3) - CN(3)*Q(2) QN(2) = CN(3)*Q(1) - CN(1)*Q(3) QN(3) = CN(1)*Q(2) - CN(2)*Q(1) ! ! Initialize parameters for the boundary traversal. ! NI = 0 EVEN = .TRUE. BP = -2. BQ = -2. PINR = .TRUE. QINR = .TRUE. I2 = LISTV(N) IF (I2 .LT. 1 .OR. I2 .GT. IMX) GO TO 12 LFT2 = CN(1)*XV(I2) + CN(2)*YV(I2) + CN(3)*ZV(I2) .GT. 0. ! ! Loop on boundary arcs I1->I2. ! DO 2 K = 1,N I1 = I2 LFT1 = LFT2 I2 = LISTV(K) IF (I2 .LT. 1 .OR. I2 .GT. IMX) GO TO 12 LFT2 = CN(1)*XV(I2) + CN(2)*YV(I2) + CN(3)*ZV(I2) .GT. 0. IF (LFT1 .EQV. LFT2) GO TO 2 ! ! I1 and I2 are on opposite sides of Q->P. Compute the ! point of intersection B. ! NI = NI + 1 V1(1) = XV(I1) V1(2) = YV(I1) V1(3) = ZV(I1) V2(1) = XV(I2) V2(2) = YV(I2) V2(3) = ZV(I2) CALL INTRSC (V1,V2,CN, B,IERR) ! ! B is between Q and P (on the shorter arc) iff ! B Forward Q->P and B Forward P->Q iff ! > 0 and > 0. ! IF (B(1)*QN(1) + B(2)*QN(2) + B(3)*QN(3) .GT. 0. & & .AND. & & B(1)*PN(1) + B(2)*PN(2) + B(3)*PN(3) .GT. 0.) & & THEN ! ! Update EVEN, BQ, QINR, BP, and PINR. ! EVEN = .NOT. EVEN D = B(1)*Q(1) + B(2)*Q(2) + B(3)*Q(3) IF (D .GT. BQ) THEN BQ = D QINR = LFT2 ENDIF D = B(1)*P(1) + B(2)*P(2) + B(3)*P(3) IF (D .GT. BP) THEN BP = D PINR = LFT1 ENDIF ENDIF 2 CONTINUE ! ! Test for consistency: NI must be even and QINR must be ! TRUE. ! IF (NI .NE. 2*(NI/2) .OR. .NOT. QINR) GO TO 1 ! ! Test for error 3: different values of PINR and EVEN. ! IF (PINR .NEQV. EVEN) GO TO 13 ! ! No error encountered. ! IER = 0 INSIDE = EVEN RETURN ! ! LV or NV is outside its valid range. ! 11 IER = 1 RETURN ! ! A LISTV entry is outside its valid range. ! 12 IER = 2 RETURN ! ! The polygon boundary is self-intersecting. ! 13 IER = 3 RETURN ! ! Consistency tests failed for all values of Q. ! 14 IER = 4 RETURN END SUBROUTINE INTADD (KK,I1,I2,I3, LIST,LPTR,LEND,LNEW ) INTEGER KK, I1, I2, I3, LIST(*), LPTR(*), LEND(*),LNEW ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/17/96 ! ! This subroutine adds an interior node to a triangulation ! of a set of points on the unit sphere. The data structure ! is updated with the insertion of node KK into the triangle ! whose vertices are I1, I2, and I3. No optimization of the ! triangulation is performed. ! ! This routine is identical to the similarly named routine ! in TRIPACK. ! ! ! On input: ! ! KK = Index of the node to be inserted. KK .GE. 1 ! and KK must not be equal to I1, I2, or I3. ! ! I1,I2,I3 = Indexes of the counterclockwise-ordered ! sequence of vertices of a triangle which ! contains node KK. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LEND,LNEW = Data structure defining the ! triangulation. Refer to Sub- ! routine TRMESH. Triangle ! (I1,I2,I3) must be included ! in the triangulation. ! ! On output: ! ! LIST,LPTR,LEND,LNEW = Data structure updated with ! the addition of node KK. KK ! will be connected to nodes I1, ! I2, and I3. ! ! Modules required by INTADD: INSERT, LSTPTR ! !*********************************************************** ! INTEGER LSTPTR INTEGER K, LP, N1, N2, N3 ! ! Local parameters: ! ! K = Local copy of KK ! LP = LIST pointer ! N1,N2,N3 = Local copies of I1, I2, and I3 ! K = KK ! ! Initialization. ! N1 = I1 N2 = I2 N3 = I3 ! ! Add K as a neighbor of I1, I2, and I3. ! LP = LSTPTR(LEND(N1),N2,LIST,LPTR) CALL INSERT (K,LP, LIST,LPTR,LNEW ) LP = LSTPTR(LEND(N2),N3,LIST,LPTR) CALL INSERT (K,LP, LIST,LPTR,LNEW ) LP = LSTPTR(LEND(N3),N1,LIST,LPTR) CALL INSERT (K,LP, LIST,LPTR,LNEW ) ! ! Add I1, I2, and I3 as neighbors of K. ! LIST(LNEW) = N1 LIST(LNEW+1) = N2 LIST(LNEW+2) = N3 LPTR(LNEW) = LNEW + 1 LPTR(LNEW+1) = LNEW + 2 LPTR(LNEW+2) = LNEW LEND(K) = LNEW + 2 LNEW = LNEW + 3 RETURN END SUBROUTINE INTRSC (P1,P2,CN, P,IER) INTEGER IER REAL P1(3), P2(3), CN(3), P(3) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/19/90 ! ! Given a great circle C and points P1 and P2 defining an ! arc A on the surface of the unit sphere, where A is the ! shorter of the two portions of the great circle C12 assoc- ! iated with P1 and P2, this subroutine returns the point ! of intersection P between C and C12 that is closer to A. ! Thus, if P1 and P2 lie in opposite hemispheres defined by ! C, P is the point of intersection of C with A. ! ! ! On input: ! ! P1,P2 = Arrays of length 3 containing the Cartesian ! coordinates of unit vectors. ! ! CN = Array of length 3 containing the Cartesian ! coordinates of a nonzero vector which defines C ! as the intersection of the plane whose normal ! is CN with the unit sphere. Thus, if C is to ! be the great circle defined by P and Q, CN ! should be P X Q. ! ! The above parameters are not altered by this routine. ! ! P = Array of length 3. ! ! On output: ! ! P = Point of intersection defined above unless IER ! .NE. 0, in which case P is not altered. ! ! IER = Error indicator. ! IER = 0 if no errors were encountered. ! IER = 1 if = . This occurs ! iff P1 = P2 or CN = 0 or there are ! two intersection points at the same ! distance from A. ! IER = 2 if P2 = -P1 and the definition of A is ! therefore ambiguous. ! ! Modules required by INTRSC: None ! ! Intrinsic function called by INTRSC: SQRT ! !*********************************************************** ! INTEGER I REAL D1, D2, PP(3), PPN, T ! ! Local parameters: ! ! D1 = ! D2 = ! I = DO-loop index ! PP = P1 + T*(P2-P1) = Parametric representation of the ! line defined by P1 and P2 ! PPN = Norm of PP ! T = D1/(D1-D2) = Parameter value chosen so that PP lies ! in the plane of C ! D1 = CN(1)*P1(1) + CN(2)*P1(2) + CN(3)*P1(3) D2 = CN(1)*P2(1) + CN(2)*P2(2) + CN(3)*P2(3) ! IF (D1 .EQ. D2) THEN IER = 1 RETURN ENDIF ! ! Solve for T such that = 0 and compute PP and PPN. ! T = D1/(D1-D2) PPN = 0. DO 1 I = 1,3 PP(I) = P1(I) + T*(P2(I)-P1(I)) PPN = PPN + PP(I)*PP(I) 1 CONTINUE ! ! PPN = 0 iff PP = 0 iff P2 = -P1 (and T = .5). ! IF (PPN .EQ. 0.) THEN IER = 2 RETURN ENDIF PPN = SQRT(PPN) ! ! Compute P = PP/PPN. ! DO 2 I = 1,3 P(I) = PP(I)/PPN 2 CONTINUE IER = 0 RETURN END INTEGER FUNCTION JRAND (N, IX,IY,IZ ) INTEGER N, IX, IY, IZ ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/28/98 ! ! This function returns a uniformly distributed pseudo- ! random integer in the range 1 to N. ! ! ! On input: ! ! N = Maximum value to be returned. ! ! N is not altered by this function. ! ! IX,IY,IZ = Integer seeds initialized to values in ! the range 1 to 30,000 before the first ! call to JRAND, and not altered between ! subsequent calls (unless a sequence of ! random numbers is to be repeated by ! reinitializing the seeds). ! ! On output: ! ! IX,IY,IZ = Updated integer seeds. ! ! JRAND = Random integer in the range 1 to N. ! ! Reference: B. A. Wichmann and I. D. Hill, "An Efficient ! and Portable Pseudo-random Number Generator", ! Applied Statistics, Vol. 31, No. 2, 1982, ! pp. 188-190. ! ! Modules required by JRAND: None ! ! Intrinsic functions called by JRAND: INT, MOD, REAL ! !*********************************************************** ! REAL U, X ! ! Local parameters: ! ! U = Pseudo-random number uniformly distributed in the ! interval (0,1). ! X = Pseudo-random number in the range 0 to 3 whose frac- ! tional part is U. ! IX = MOD(171*IX,30269) IY = MOD(172*IY,30307) IZ = MOD(170*IZ,30323) X = (REAL(IX)/30269.) + (REAL(IY)/30307.) + (REAL(IZ)/30323.) U = X - INT(X) JRAND = REAL(N)*U + 1. RETURN END LOGICAL FUNCTION LEFT (X1,Y1,Z1,X2,Y2,Z2,X0,Y0,Z0) REAL X1, Y1, Z1, X2, Y2, Z2, X0, Y0, Z0 ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/15/96 ! ! This function determines whether node N0 is in the ! (closed) left hemisphere defined by the plane containing ! N1, N2, and the origin, where left is defined relative to ! an observer at N1 facing N2. ! ! ! On input: ! ! X1,Y1,Z1 = Coordinates of N1. ! ! X2,Y2,Z2 = Coordinates of N2. ! ! X0,Y0,Z0 = Coordinates of N0. ! ! Input parameters are not altered by this function. ! ! On output: ! ! LEFT = TRUE if and only if N0 is in the closed ! left hemisphere. ! ! Modules required by LEFT: None ! !*********************************************************** ! ! LEFT = TRUE iff = det(N0,N1,N2) .GE. 0. ! LEFT = X0*(Y1*Z2-Y2*Z1) - Y0*(X1*Z2-X2*Z1) + & & Z0*(X1*Y2-X2*Y1) .GE. 0. RETURN END INTEGER FUNCTION LSTPTR (LPL,NB,LIST,LPTR) INTEGER LPL, NB, LIST(*), LPTR(*) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/15/96 ! ! 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. ! ! ! On input: ! ! LPL = LEND(N0) ! ! NB = Index of the node whose pointer is to be re- ! turned. NB must be connected to N0. ! ! LIST,LPTR = Data structure defining the triangula- ! tion. Refer to Subroutine TRMESH. ! ! Input parameters are not altered by this function. ! ! On output: ! ! LSTPTR = Pointer such that LIST(LSTPTR) = NB or ! LIST(LSTPTR) = -NB, unless NB is not a ! neighbor of N0, in which case LSTPTR = LPL. ! ! Modules required by LSTPTR: None ! !*********************************************************** ! INTEGER LP, ND ! ! Local parameters: ! ! LP = LIST pointer ! ND = Nodal index ! LP = LPTR(LPL) 1 ND = LIST(LP) IF (ND .EQ. NB) GO TO 2 LP = LPTR(LP) IF (LP .NE. LPL) GO TO 1 ! 2 LSTPTR = LP RETURN END INTEGER FUNCTION NBCNT (LPL,LPTR) INTEGER LPL, LPTR(*) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/15/96 ! ! This function returns the number of neighbors of a node ! N0 in a triangulation created by Subroutine TRMESH. ! ! This function is identical to the similarly named ! function in TRIPACK. ! ! ! On input: ! ! LPL = LIST pointer to the last neighbor of N0 -- ! LPL = LEND(N0). ! ! LPTR = Array of pointers associated with LIST. ! ! Input parameters are not altered by this function. ! ! On output: ! ! NBCNT = Number of neighbors of N0. ! ! Modules required by NBCNT: None ! !*********************************************************** ! INTEGER K, LP ! ! Local parameters: ! ! K = Counter for computing the number of neighbors ! LP = LIST pointer ! LP = LPL K = 1 ! 1 LP = LPTR(LP) IF (LP .EQ. LPL) GO TO 2 K = K + 1 GO TO 1 ! 2 NBCNT = K RETURN END INTEGER FUNCTION NEARND (P,IST,N,X,Y,Z,LIST,LPTR, LEND, AL) INTEGER IST, N, LIST(*), LPTR(*), LEND(N) REAL P(3), X(N), Y(N), Z(N), AL ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/28/98 ! ! Given a point P on the surface of the unit sphere and a ! Delaunay triangulation created by Subroutine 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. ! ! ! On input: ! ! P = Array of length 3 containing the Cartesian coor- ! dinates of the point P to be located relative to ! the triangulation. It is assumed without a test ! that P(1)**2 + P(2)**2 + P(3)**2 = 1. ! ! IST = Index of a node at which TRFIND begins the ! search. Search time depends on the proximity ! of this node to P. ! ! N = Number of nodes in the triangulation. N .GE. 3. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes. ! ! LIST,LPTR,LEND = Data structure defining the trian- ! gulation. Refer to TRMESH. ! ! Input parameters are not altered by this function. ! ! On output: ! ! NEARND = Nodal index of the nearest node to P, or 0 ! if N < 3 or the triangulation data struc- ! ture is invalid. ! ! AL = Arc length (angular distance in radians) be- ! tween P and NEARND unless NEARND = 0. ! ! Note that the number of candidates for NEARND ! (neighbors of P) is limited to LMAX defined in ! the PARAMETER statement below. ! ! Modules required by NEARND: JRAND, LSTPTR, TRFIND, STORE ! ! Intrinsic functions called by NEARND: ABS, ACOS ! !*********************************************************** ! INTEGER LSTPTR INTEGER LMAX PARAMETER (LMAX=25) INTEGER I1, I2, I3, L, LISTP(LMAX), LP, LP1, LP2, & & LPL, LPTRP(LMAX), N1, N2, N3, NN, NR, NST REAL B1, B2, B3, DS1, DSR, DX1, DX2, DX3, DY1, & & DY2, DY3, DZ1, DZ2, DZ3 ! ! 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 ! ! ! Store local parameters and test for N invalid. ! NN = N IF (NN .LT. 3) GO TO 6 NST = IST IF (NST .LT. 1 .OR. NST .GT. NN) NST = 1 ! ! 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 .EQ. 0) GO TO 6 ! ! 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 .NE. 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. ! 1 LPL = LEND(N1) N1 = -LIST(LPL) L = LP1 LP1 = L+1 LISTP(L) = N1 LPTRP(L) = LP1 IF (N1 .NE. I2 .AND. LP1 .LT. LMAX) GO TO 1 L = LP1 LISTP(L) = 0 LPTRP(L) = 1 ENDIF ! ! 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. ! 2 LP = LSTPTR(LEND(N1),N2,LIST,LPTR) IF (LIST(LP) .LT. 0) GO TO 3 LP = LPTR(LP) N3 = ABS(LIST(LP)) ! ! Swap test: Exit the loop if L = LMAX. ! IF (L .EQ. LMAX) GO TO 4 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) IF ( DX3*(DY2*DZ1 - DY1*DZ2) - & & DY3*(DX2*DZ1 - DX1*DZ2) + & & DZ3*(DX2*DY1 - DX1*DY2) .LE. 0. ) GO TO 3 ! ! Swap: Insert N3 following N2 in the adjacency list for P. ! The two new arcs opposite P must be tested. ! L = L+1 LPTRP(LP2) = L LISTP(L) = N3 LPTRP(L) = LP1 LP1 = L N1 = N3 GO TO 2 ! ! No swap: Advance to the next arc and test for termination ! on N1 = I1 (LP1 = 1) or N1 followed by 0. ! 3 IF (LP1 .EQ. 1) GO TO 4 LP2 = LP1 N2 = N1 LP1 = LPTRP(LP1) N1 = LISTP(LP1) IF (N1 .EQ. 0) GO TO 4 GO TO 2 ! ! 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. ! 4 NR = I1 DSR = -(X(NR)*P(1) + Y(NR)*P(2) + Z(NR)*P(3)) DO 5 LP = 2,L N1 = LISTP(LP) IF (N1 .EQ. 0) GO TO 5 DS1 = -(X(N1)*P(1) + Y(N1)*P(2) + Z(N1)*P(3)) IF (DS1 .LT. DSR) THEN NR = N1 DSR = DS1 ENDIF 5 CONTINUE DSR = -DSR IF (DSR .GT. 1.0) DSR = 1.0 AL = ACOS(DSR) NEARND = NR RETURN ! ! Invalid input. ! 6 NEARND = 0 RETURN END SUBROUTINE OPTIM (X,Y,Z,NA, LIST,LPTR,LEND,NIT, IWK, IER) INTEGER NA, LIST(*), LPTR(*), LEND(*), NIT, IWK(2,NA), IER REAL X(*), Y(*), Z(*) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/30/98 ! ! 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. ! ! ! On input: ! ! X,Y,Z = Arrays containing the nodal coordinates. ! ! NA = Number of arcs in the set. NA .GE. 0. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LEND = Data structure defining the trian- ! gulation. Refer to Subroutine ! TRMESH. ! ! NIT = Maximum number of iterations to be performed. ! NIT = 4*NA should be sufficient. NIT .GE. 1. ! ! IWK = Integer array dimensioned 2 by NA containing ! the nodal indexes of the arc endpoints (pairs ! of endpoints are stored in columns). ! ! On output: ! ! LIST,LPTR,LEND = Updated triangulation data struc- ! ture reflecting the swaps. ! ! NIT = Number of iterations performed. ! ! IWK = Endpoint indexes of the new set of arcs ! reflecting the swaps. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 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. ! IER = 2 if NA < 0 or NIT < 1 on input. ! IER = 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. ! IER = 4 if a zero pointer was returned by ! Subroutine SWAP. ! ! Modules required by OPTIM: LSTPTR, SWAP, SWPTST ! ! Intrinsic function called by OPTIM: ABS ! !*********************************************************** ! INTEGER I, IO1, IO2, ITER, LP, LP21, LPL, LPP, MAXIT, N1, N2, NNA LOGICAL SWPTST LOGICAL SWP ! ! 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 ! NNA = NA MAXIT = NIT IF (NNA .LT. 0 .OR. MAXIT .LT. 1) GO TO 7 ! ! Initialize iteration count ITER and test for NA = 0. ! ITER = 0 IF (NNA .EQ. 0) GO TO 5 ! ! Top of loop -- ! SWP = TRUE iff a swap occurred in the current iteration. ! 1 IF (ITER .EQ. MAXIT) GO TO 6 ITER = ITER + 1 SWP = .FALSE. ! ! Inner loop on arcs IO1-IO2 -- ! DO 4 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) 2 IF (LIST(LP) .EQ. IO2) GO TO 3 LPP = LP LP = LPTR(LPP) IF (LP .NE. LPL) GO TO 2 ! ! 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)) .NE. IO2) GO TO 8 IF (LIST(LP) .LT. 0) GO TO 4 ! ! Store N1 and N2, or bypass the swap test if IO1 is a ! boundary node and IO2 is its first neighbor. ! 3 N2 = LIST(LPP) IF (N2 .LT. 0) GO TO 4 LP = LPTR(LP) N1 = ABS(LIST(LP)) ! ! Test IO1-IO2 for a swap, and update IWK if necessary. ! IF ( .NOT. SWPTST(N1,N2,IO1,IO2,X,Y,Z) ) GO TO 4 CALL SWAP (N1,N2,IO1,IO2, LIST,LPTR,LEND, LP21) IF (LP21 .EQ. 0) GO TO 9 SWP = .TRUE. IWK(1,I) = N1 IWK(2,I) = N2 4 CONTINUE IF (SWP) GO TO 1 ! ! Successful termination. ! 5 NIT = ITER IER = 0 RETURN ! ! MAXIT iterations performed without convergence. ! 6 NIT = MAXIT IER = 1 RETURN ! ! Invalid input parameter. ! 7 NIT = 0 IER = 2 RETURN ! ! IO2 is not a neighbor of IO1. ! 8 NIT = ITER IER = 3 RETURN ! ! Zero pointer returned by SWAP. ! 9 NIT = ITER IER = 4 RETURN END SUBROUTINE SCOORD (PX,PY,PZ, PLAT,PLON,PNRM) REAL PX, PY, PZ, PLAT, PLON, PNRM ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 08/27/90 ! ! This subroutine converts a point P from Cartesian coor- ! dinates to spherical coordinates. ! ! ! On input: ! ! PX,PY,PZ = Cartesian coordinates of P. ! ! Input parameters are not altered by this routine. ! ! On output: ! ! PLAT = Latitude of P in the range -PI/2 to PI/2, or ! 0 if PNRM = 0. PLAT should be scaled by ! 180/PI to obtain the value in degrees. ! ! PLON = Longitude of P in the range -PI to PI, or 0 ! if P lies on the Z-axis. PLON should be ! scaled by 180/PI to obtain the value in ! degrees. ! ! PNRM = Magnitude (Euclidean norm) of P. ! ! Modules required by SCOORD: None ! ! Intrinsic functions called by SCOORD: ASIN, ATAN2, SQRT ! !*********************************************************** ! PNRM = SQRT(PX*PX + PY*PY + PZ*PZ) IF (PX .NE. 0. .OR. PY .NE. 0.) THEN PLON = ATAN2(PY,PX) ELSE PLON = 0. ENDIF IF (PNRM .NE. 0.) THEN PLAT = ASIN(PZ/PNRM) ELSE PLAT = 0. ENDIF RETURN END REAL FUNCTION STORE (X) REAL X ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 05/09/92 ! ! 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. ! ! ! On input: ! ! X = Value to be stored. ! ! X is not altered by this function. ! ! On output: ! ! STORE = Value of X after it has been stored and ! possibly truncated or rounded to the single ! precision word length. ! ! Modules required by STORE: None ! !*********************************************************** ! REAL Y COMMON/STCOM/Y Y = X STORE = Y RETURN END SUBROUTINE SWAP (IN1,IN2,IO1,IO2, LIST,LPTR, LEND, LP21) INTEGER IN1, IN2, IO1, IO2, LIST(*), LPTR(*), LEND(*), LP21 ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 06/22/98 ! ! 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. ! ! ! On input: ! ! IN1,IN2,IO1,IO2 = Nodal indexes of the vertices of ! the quadrilateral. IO1-IO2 is re- ! placed by IN1-IN2. (IO1,IO2,IN1) ! and (IO2,IO1,IN2) must be trian- ! gles on input. ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR,LEND = Data structure defining the trian- ! gulation. Refer to Subroutine ! TRMESH. ! ! On output: ! ! LIST,LPTR,LEND = Data structure updated with the ! swap -- triangles (IO1,IO2,IN1) and ! (IO2,IO1,IN2) are replaced by ! (IN1,IN2,IO2) and (IN2,IN1,IO1) ! unless LP21 = 0. ! ! 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. ! ! Module required by SWAP: LSTPTR ! ! Intrinsic function called by SWAP: ABS ! !*********************************************************** ! INTEGER LSTPTR INTEGER LP, LPH, LPSAV ! ! Local parameters: ! ! LP,LPH,LPSAV = LIST pointers ! ! ! Test for IN1 and IN2 adjacent. ! LP = LSTPTR(LEND(IN1),IN2,LIST,LPTR) IF (ABS(LIST(LP)) .EQ. IN2) THEN LP21 = 0 RETURN ENDIF ! ! 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) .EQ. LPH) LEND(IO1) = LP ! ! 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) .EQ. LPH) LEND(IO2) = LP ! ! 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 LOGICAL FUNCTION SWPTST (N1,N2,N3,N4,X,Y,Z) INTEGER N1, N2, N3, N4 REAL X(*), Y(*), Z(*) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 03/29/91 ! ! 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. ! ! ! On input: ! ! 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 connec- ! ting N3 to N4 if SWPTST = TRUE. Refer ! to Subroutine SWAP. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes. (X(I),Y(I),Z(I)) ! define node I for I = N1, N2, N3, and N4. ! ! Input parameters are not altered by this routine. ! ! On output: ! ! SWPTST = TRUE if and only if the arc connecting N1 ! and N2 should be swapped for an arc con- ! necting N3 and N4. ! ! Modules required by SWPTST: None ! !*********************************************************** ! REAL DX1, DX2, DX3, DY1, DY2, DY3, DZ1, DZ2, DZ3, X4, Y4, Z4 ! ! 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 ! 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) .GT. 0. RETURN END SUBROUTINE TRANS (N,RLAT,RLON, X,Y,Z) INTEGER N REAL RLAT(N), RLON(N), X(N), Y(N), Z(N) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 04/08/90 ! ! This subroutine transforms spherical coordinates into ! Cartesian coordinates on the unit sphere for input to ! Subroutine TRMESH. Storage for X and Y may coincide with ! storage for RLAT and RLON if the latter need not be saved. ! ! ! On input: ! ! N = Number of nodes (points on the unit sphere) ! whose coordinates are to be transformed. ! ! RLAT = Array of length N containing latitudinal ! coordinates of the nodes in radians. ! ! RLON = Array of length N containing longitudinal ! coordinates of the nodes in radians. ! ! The above parameters are not altered by this routine. ! ! X,Y,Z = Arrays of length at least N. ! ! On output: ! ! X,Y,Z = Cartesian coordinates in the range -1 to 1. ! X(I)**2 + Y(I)**2 + Z(I)**2 = 1 for I = 1 ! to N. ! ! Modules required by TRANS: None ! ! Intrinsic functions called by TRANS: COS, SIN ! !*********************************************************** ! INTEGER I, NN REAL COSPHI, PHI, THETA ! ! Local parameters: ! ! COSPHI = cos(PHI) ! I = DO-loop index ! NN = Local copy of N ! PHI = Latitude ! THETA = Longitude ! NN = N DO 1 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) 1 CONTINUE RETURN END SUBROUTINE TRFIND(NST,P,N,X,Y,Z,LIST,LPTR,LEND, B1,B2,B3,I1,I2,I3) INTEGER NST, N, LIST(*), LPTR(*), LEND(N), I1, I2, I3 REAL P(3), X(N), Y(N), Z(N), B1, B2, B3 ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 11/30/99 ! ! This subroutine locates a point P relative to a triangu- ! lation created by Subroutine TRMESH. If P is contained in ! a triangle, the three vertex indexes and barycentric coor- ! dinates are returned. Otherwise, the indexes of the ! visible boundary nodes are returned. ! ! ! On input: ! ! NST = Index of a node at which TRFIND begins its ! search. Search time depends on the proximity ! of this node to P. ! ! P = Array of length 3 containing the x, y, and z ! coordinates (in that order) of the point P to be ! located. ! ! N = Number of nodes in the triangulation. N .GE. 3. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the triangulation nodes (unit ! vectors). (X(I),Y(I),Z(I)) defines node I ! for I = 1 to N. ! ! LIST,LPTR,LEND = Data structure defining the trian- ! gulation. Refer to Subroutine ! TRMESH. ! ! Input parameters are not altered by this routine. ! ! On output: ! ! B1,B2,B3 = Unnormalized barycentric coordinates of ! the central projection of P onto the un- ! derlying planar triangle if P is in the ! convex hull of the nodes. These parame- ! ters are not altered if I1 = 0. ! ! I1,I2,I3 = Counterclockwise-ordered vertex indexes ! of a triangle containing P if P is con- ! tained 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 vis- ! ible 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. ! ! Modules required by TRFIND: JRAND, LSTPTR, STORE ! ! Intrinsic function called by TRFIND: ABS ! !*********************************************************** ! INTEGER JRAND, LSTPTR INTEGER IX, IY, IZ, LP, N0, N1, N1S, N2, N2S, N3, N4, NEXT, NF, NL REAL STORE REAL DET, EPS, PTN1, PTN2, Q(3), S12, TOL, XP, YP, ZP REAL X0, X1, X2, Y0, Y1, Y2, Z0, Z1, Z2 REAL DI1PSQ, DI2PSQ, DI1I2SQ ! SAVE IX, IY, IZ DATA IX/1/, IY/2/, IZ/3/ ! ! 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 ! ! DI1PSQ, DI2PSQ, DI1I2SQ - from point outside triangulation to ! two boundary points. Used in triangulation ! for finding boundary weighted interpolation ! ! Statement function: ! ! DET(X1,...,Z0) .GE. 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 ob- ! server 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 .LT. 1 .OR. N0 .GT. N) & & N0 = JRAND(N, IX,IY,IZ ) ! ! Compute the relative machine precision EPS and TOL. ! EPS = 1.E0 1 EPS = EPS/2.E0 IF (STORE(EPS+1.E0) .GT. 1.E0) GO TO 1 EPS = 2.E0*EPS TOL = 100.E0*EPS ! ! Set NF and NL to the first and last neighbors of N0, and ! initialize N1 = NF. ! 2 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 (NL .GT. 0) THEN ! ! N0 is an interior node. Find N1. ! 3 IF ( DET(X(N0),Y(N0),Z(N0),X(N1),Y(N1),Z(N1), & & XP,YP,ZP) .LT. 0. ) THEN LP = LPTR(LP) N1 = LIST(LP) IF (N1 .EQ. NL) GO TO 6 GO TO 3 ENDIF ELSE ! ! N0 is a boundary node. Test for P exterior. ! NL = -NL IF ( DET(X(N0),Y(N0),Z(N0),X(NF),Y(NF),Z(NF), & & XP,YP,ZP) .LT. 0. ) THEN ! ! P is to the right of the boundary edge N0->NF. ! N1 = N0 N2 = NF GO TO 9 ENDIF IF ( DET(X(NL),Y(NL),Z(NL),X(N0),Y(N0),Z(N0), & & XP,YP,ZP) .LT. 0. ) THEN ! ! P is to the right of the boundary edge NL->N0. ! N1 = NL N2 = N0 GO TO 9 ENDIF ENDIF ! ! P is to the left of arcs N0->N1 and NL->N0. Set N2 to the ! next neighbor of N0 (following N1). ! 4 LP = LPTR(LP) N2 = ABS(LIST(LP)) IF ( DET(X(N0),Y(N0),Z(N0),X(N2),Y(N2),Z(N2), & & XP,YP,ZP) .LT. 0. ) GO TO 7 N1 = N2 IF (N1 .NE. NL) GO TO 4 IF ( DET(X(N0),Y(N0),Z(N0),X(NF),Y(NF),Z(NF), & & XP,YP,ZP) .LT. 0. ) GO TO 6 ! ! 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)) & & .LT. 1.0-4.0*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. ! 5 IF ( DET(X(N1),Y(N1),Z(N1),X(N0),Y(N0),Z(N0), & & XP,YP,ZP) .GE. 0. ) THEN LP = LPTR(LP) N1 = ABS(LIST(LP)) IF (N1 .EQ. NL) GO TO 14 GO TO 5 ENDIF ENDIF ! ! 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 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 N3 = N0 N1S = N1 N2S = N2 ! ! Top of edge-hopping loop: ! 8 B3 = DET(X(N1),Y(N1),Z(N1),X(N2),Y(N2),Z(N2),XP,YP,ZP) IF (B3 .LT. 0.) 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) .LT. 0) GO TO 9 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) .LT. 0. ) THEN N3 = N2 N2 = N4 N1S = N1 IF (N2 .NE. N2S .AND. N2 .NE. N0) GO TO 8 ELSE N3 = N1 N1 = N4 N2S = N2 IF (N1 .NE. N1S .AND. N1 .NE. N0) GO TO 8 ENDIF ! ! 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 ENDIF ! ! P is in (N1,N2,N3) unless N0, N1, N2, and P are collinear ! or P is close to -N0. ! IF (B3 .GE. EPS) THEN ! ! B3 .NE. 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) IF (B1 .LT. -TOL .OR. B2 .LT. -TOL) THEN ! ! Restart with N0 randomly selected. ! N0 = JRAND(N, IX,IY,IZ ) GO TO 2 ENDIF 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. 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 IF (B1 .LT. -TOL .OR. B2 .LT. -TOL) THEN ! ! Restart with N0 randomly selected. ! N0 = JRAND(N, IX,IY,IZ ) GO TO 2 ENDIF ENDIF ! ! P is in (N1,N2,N3). ! I1 = N1 I2 = N2 I3 = N3 IF (B1 .LT. 0.0) B1 = 0.0 IF (B2 .LT. 0.0) B2 = 0.0 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 N1S = N1 N2S = N2 NL = 0 ! ! Counterclockwise Boundary Traversal: ! 10 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) .GE. 0. ) 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) .GE. 0.) GO TO 11 IF (X(NEXT)*Q(1) + Y(NEXT)*Q(2) + Z(NEXT)*Q(3) & & .GE. 0.) GO TO 11 ! ! N1, N2, NEXT, and P are nearly collinear, and N2 is ! the leftmost visible node. ! NL = N2 ENDIF ! ! Bottom of counterclockwise loop: ! N1 = N2 N2 = NEXT IF (N2 .NE. N1S) GO TO 10 ! ! All boundary nodes are visible from P. ! I1 = N1S I2 = N1S I3 = 0 RETURN ! ! N2 is the rightmost visible node. ! 11 NF = N2 IF (NL .EQ. 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 LP = LEND(N1) NEXT = -LIST(LP) IF ( DET(X(NEXT),Y(NEXT),Z(NEXT),X(N1),Y(N1),Z(N1), & & XP,YP,ZP) .GE. 0. ) 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) .GE. 0.) GO TO 13 IF (X(NEXT)*Q(1) + Y(NEXT)*Q(2) + Z(NEXT)*Q(3) & & .GE. 0.) GO TO 13 ! ! P, NEXT, N1, and N2 are nearly collinear and N1 is the ! rightmost visible node. ! NF = N1 ENDIF ! ! Bottom of clockwise loop: ! N2 = N1 N1 = NEXT IF (N1 .NE. N1S) GO TO 12 ! ! All boundary nodes are visible from P. ! I1 = N1 I2 = N1 I3 = 0 ! added to give valid return for two points B1 = 0.5 B2 = 0.5 RETURN ! ! N1 is the leftmost visible node. ! 13 NL = N1 ENDIF ! ! NF and NL have been found. ! I1 = NF I2 = NL I3 = 0 ! added to give valid return for two points ! this is a linear interpolation between the two points using ! the perpendiculars to the line between the two points as the ! limit. Outside this limit, the end point is given full weight ! Length of the three sides (2 boundary points and 1 trial point) ! and law of cosines used to find distance of trial point between ! the two boundary points ! Length of three sides DI1PSQ = (P(1)-X(I1))**2 + (P(2)-Y(I1))**2 + (P(3)-Z(I1))**2 DI2PSQ = (P(1)-X(I2))**2 + (P(2)-Y(I2))**2 + (P(3)-Z(I2))**2 DI1I2SQ = (X(I1)-X(I2))**2 + (Y(I1)-Y(I2))**2 + (Z(I1)-Z(I2))**2 ! weighted distance of trial point between two boundary points B1 = (DI2PSQ + DI1I2SQ - DI1PSQ) / (2*DI1I2SQ) B2 = (DI1PSQ + DI1I2SQ - DI2PSQ) / (2*DI1I2SQ) IF( B1 .LT. 0.0 ) THEN ! point is outside perpendicular rectangle B1 = 0.0 B2 = 1.0 ELSE IF( B2 .LT. 0.0 ) THEN ! point is outside perpendicular rectangle B1 = 1.0 B2 = 0.0 END IF RETURN ! ! All points are collinear (coplanar). ! 14 I1 = 0 I2 = 0 I3 = 0 RETURN END SUBROUTINE TRLIST (N,LIST,LPTR,LEND,NROW, NT,LTRI,IER) INTEGER N, LIST(*), LPTR(*), LEND(N), NROW, NT, LTRI(NROW,*), IER ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/20/96 ! ! This subroutine converts a triangulation data structure ! from the linked list created by Subroutine TRMESH to a ! triangle list. ! ! On input: ! ! N = Number of nodes in the triangulation. N .GE. 3. ! ! LIST,LPTR,LEND = Linked list data structure defin- ! ing the triangulation. Refer to ! Subroutine TRMESH. ! ! NROW = Number of rows (entries per triangle) re- ! served 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. ! ! The above parameters are not altered by this routine. ! ! LTRI = Integer array of length at least NROW*NT, ! where NT is at most 2N-4. (A sufficient ! length is 12N if NROW=6 or 18N if NROW=9.) ! ! On output: ! ! NT = Number of triangles in the triangulation unless ! IER .NE. 0, in which case NT = 0. NT = 2N-NB-2 ! if NB .GE. 3 or 2N-4 if NB = 0, where NB is the ! number of boundary nodes. ! ! LTRI = NROW by NT array whose 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 .GE. 3 ! or 3N-6 if NB = 0. The triangles are or- ! dered on first (smallest) vertex indexes. ! ! IER = Error indicator. ! IER = 0 if no errors were encountered. ! IER = 1 if N or NROW is outside its valid ! range on input. ! IER = 2 if the triangulation data structure ! (LIST,LPTR,LEND) is invalid. Note, ! however, that these arrays are not ! completely tested for validity. ! ! Modules required by TRLIST: None ! ! Intrinsic function called by TRLIST: ABS ! !*********************************************************** ! INTEGER I, I1, I2, I3, ISV, J, KA, KN, KT, LP, LP2, & & LPL, LPLN1, N1, N2, N3, NM2 LOGICAL ARCS ! ! 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 ! ! ! Test for invalid input parameters. ! IF (N .LT. 3 .OR. (NROW .NE. 6 .AND. NROW .NE. 9)) & & GO TO 11 ! ! 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 .EQ. 9 KA = 0 KT = 0 NM2 = N-2 ! ! Loop on nodes N1. ! DO 9 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 LP2 = LPTR(LP2) N2 = LIST(LP2) LP = LPTR(LP2) N3 = ABS(LIST(LP)) IF (N2 .LT. N1 .OR. N3 .LT. N1) GO TO 8 ! ! 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 7 I = 1,3 IF (I .EQ. 1) THEN I1 = N3 I2 = N2 ELSEIF (I .EQ. 2) THEN I1 = N1 I2 = N3 ELSE I1 = N2 I2 = N1 ENDIF ! ! Set I3 to the neighbor of I1 that follows I2 unless ! I2->I1 is a boundary arc. ! LPL = LEND(I1) LP = LPTR(LPL) 2 IF (LIST(LP) .EQ. I2) GO TO 3 LP = LPTR(LP) IF (LP .NE. LPL) GO TO 2 ! ! I2 is the last neighbor of I1 unless the data structure ! is invalid. Bypass the search for a neighboring ! triangle if I2->I1 is a boundary arc. ! IF (ABS(LIST(LP)) .NE. I2) GO TO 12 KN = 0 IF (LIST(LP) .LT. 0) GO TO 6 ! ! I2->I1 is not a boundary arc, and LP points to I2 as ! a neighbor of I1. ! 3 LP = LPTR(LP) I3 = ABS(LIST(LP)) ! ! Find J such that LTRI(J,KN) = I3 (not used if KN > KT), ! and permute the vertex indexes of KN so that I1 is ! smallest. ! IF (I1 .LT. I2 .AND. I1 .LT. I3) THEN J = 3 ELSEIF (I2 .LT. I3) THEN J = 2 ISV = I1 I1 = I2 I2 = I3 I3 = ISV ELSE J = 1 ISV = I1 I1 = I3 I3 = I2 I2 = ISV ENDIF ! ! Test for KN > KT (triangle index not yet assigned). ! IF (I1 .GT. N1) GO TO 7 ! ! Find KN, if it exists, by searching the triangle list in ! reverse order. ! DO 4 KN = KT-1,1,-1 IF (LTRI(1,KN) .EQ. I1 .AND. LTRI(2,KN) .EQ. & & I2 .AND. LTRI(3,KN) .EQ. I3) GO TO 5 4 CONTINUE GO TO 7 ! ! Store KT as a neighbor of KN. ! 5 LTRI(J+3,KN) = KT ! ! Store KN as a neighbor of KT, and add a new arc KA. ! 6 LTRI(I+3,KT) = KN IF (ARCS) THEN KA = KA + 1 LTRI(I+6,KT) = KA IF (KN .NE. 0) LTRI(J+6,KN) = KA ENDIF 7 CONTINUE ! ! Bottom of loop on triangles. ! 8 IF (LP2 .NE. LPLN1) GO TO 1 9 CONTINUE ! ! No errors encountered. ! NT = KT IER = 0 RETURN ! ! Invalid input parameter. ! 11 NT = 0 IER = 1 RETURN ! ! Invalid triangulation data structure: I1 is a neighbor of ! I2, but I2 is not a neighbor of I1. ! 12 NT = 0 IER = 2 RETURN END SUBROUTINE TRLPRT (N,X,Y,Z,IFLAG,NROW,NT,LTRI,LOUT) INTEGER N, IFLAG, NROW, NT, LTRI(NROW,NT), LOUT REAL X(N), Y(N), Z(N) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/02/98 ! ! This subroutine prints the triangle list created by Sub- ! routine TRLIST and, optionally, the nodal coordinates ! (either latitude and longitude or Cartesian coordinates) ! on logical unit LOUT. The numbers of boundary nodes, ! triangles, and arcs are also printed. ! ! ! On input: ! ! N = Number of nodes in the triangulation. ! 3 .LE. N .LE. 9999. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes if IFLAG = 0, or ! (X and Y only) arrays of length N containing ! longitude and latitude, respectively, if ! IFLAG > 0, or unused dummy parameters if ! IFLAG < 0. ! ! IFLAG = Nodal coordinate option indicator: ! IFLAG = 0 if X, Y, and Z (assumed to contain ! Cartesian coordinates) are to be ! printed (to 6 decimal places). ! IFLAG > 0 if only X and Y (assumed to con- ! tain longitude and latitude) are ! to be printed (to 6 decimal ! places). ! IFLAG < 0 if only the adjacency lists are to ! be printed. ! ! NROW = Number of rows (entries per triangle) re- ! served 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. ! ! NT = Number of triangles in the triangulation. ! 1 .LE. NT .LE. 9999. ! ! LTRI = NROW by NT array whose 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. ! ! LOUT = Logical unit number for output. If LOUT is ! not in the range 0 to 99, output is written ! to unit 6. ! ! Input parameters are not altered by this routine. ! ! On output: ! ! The triangle list and nodal coordinates (as specified by ! IFLAG) are written to unit LOUT. ! ! Modules required by TRLPRT: None ! !*********************************************************** ! INTEGER I, K, LUN, NA, NB, NL, NLMAX, NMAX DATA NMAX/9999/, NLMAX/58/ ! ! Local parameters: ! ! I = DO-loop, nodal index, and row index for LTRI ! K = DO-loop and triangle index ! LUN = Logical unit number for output ! 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 addi- ! tional lines) ! NMAX = Maximum value of N and NT (4-digit format) ! LUN = LOUT IF (LUN .LT. 0 .OR. LUN .GT. 99) LUN = 6 ! ! Print a heading and test for invalid input. ! WRITE (LUN,100) N NL = 3 IF (N .LT. 3 .OR. N .GT. NMAX .OR. & & (NROW .NE. 6 .AND. NROW .NE. 9) .OR. & & NT .LT. 1 .OR. NT .GT. NMAX) THEN ! ! Print an error message and exit. ! WRITE (LUN,110) N, NROW, NT RETURN ENDIF IF (IFLAG .EQ. 0) THEN ! ! Print X, Y, and Z. ! WRITE (LUN,101) NL = 6 DO 1 I = 1,N IF (NL .GE. NLMAX) THEN WRITE (LUN,108) NL = 0 ENDIF WRITE (LUN,103) I, X(I), Y(I), Z(I) NL = NL + 1 1 CONTINUE ELSEIF (IFLAG .GT. 0) THEN ! ! Print X (longitude) and Y (latitude). ! WRITE (LUN,102) NL = 6 DO 2 I = 1,N IF (NL .GE. NLMAX) THEN WRITE (LUN,108) NL = 0 ENDIF WRITE (LUN,104) I, X(I), Y(I) NL = NL + 1 2 CONTINUE ENDIF ! ! Print the triangulation LTRI. ! IF (NL .GT. NLMAX/2) THEN WRITE (LUN,108) NL = 0 ENDIF IF (NROW .EQ. 6) THEN WRITE (LUN,105) ELSE WRITE (LUN,106) ENDIF NL = NL + 5 DO 3 K = 1,NT IF (NL .GE. NLMAX) THEN WRITE (LUN,108) NL = 0 ENDIF WRITE (LUN,107) K, (LTRI(I,K), I = 1,NROW) NL = NL + 1 3 CONTINUE ! ! Print NB, NA, and NT (boundary nodes, arcs, and ! triangles). ! NB = 2*N - NT - 2 IF (NB .LT. 3) THEN NB = 0 NA = 3*N - 6 ELSE NA = NT + N - 1 ENDIF WRITE (LUN,109) NB, NA, 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',8X,'Longitude',9X,'Latitude'//) 103 FORMAT (8X,I4,3E17.6) 104 FORMAT (16X,I4,2E17.6) 105 FORMAT (//1X,'Triangle',8X,'Vertices',12X,'Neighbors'/ & & 4X,'KT',7X,'N1',5X,'N2',5X,'N3',4X,'KT1',4X, & & 'KT2',4X,'KT3'/) 106 FORMAT (//1X,'Triangle',8X,'Vertices',12X,'Neighbors', & & 14X,'Arcs'/ & & 4X,'KT',7X,'N1',5X,'N2',5X,'N3',4X,'KT1',4X, & & 'KT2',4X,'KT3',4X,'KA1',4X,'KA2',4X,'KA3'/) 107 FORMAT (2X,I4,2X,6(3X,I4),3(2X,I5)) 108 FORMAT (///) 109 FORMAT (/1X,'NB = ',I4,' Boundary Nodes',5X, & & 'NA = ',I5,' Arcs',5X,'NT = ',I5, & & ' Triangles') 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) INTEGER N, LIST(*), LPTR(*), LEND(N), LNEW, NEAR(N), NEXT(N), IER REAL X(N), Y(N), Z(N), DIST(N) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/08/99 ! ! 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 en- ! tire 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 trian- ! gle. ! ! 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 connec- ! ted by an arc in the triangulation. ! ! GETNP - Determines the ordered sequence of L closest ! nodes to a given node, along with the associ- ! ated 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 fin- ! ite element code. ! ! TRLPRT - Prints the triangle list created by Subroutine ! 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. ! ! ! On input: ! ! N = Number of nodes in the triangulation. N .GE. 3. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of distinct nodes. (X(K),Y(K), ! Z(K)) is referred to as node K, and K is re- ! ferred 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 col- ! linear (lie on a common great circle). ! ! The above parameters are not altered by this routine. ! ! LIST,LPTR = Arrays of length at least 6N-12. ! ! LEND = Array of length at least N. ! ! NEAR,NEXT,DIST = Work space arrays of length at ! least N. The space is used to ! efficiently determine the nearest ! triangulation node to each un- ! processed node for use by ADDNOD. ! ! On output: ! ! LIST = Set of 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 bound- ! ary node are boundary nodes (the first neigh- ! bor 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. ! ! LPTR = 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 neigh- ! bor). ! ! LEND = Set of pointers to adjacency lists. LEND(K) ! points to the last neighbor of node K for ! K = 1,...,N. Thus, LIST(LEND(K)) < 0 if and ! only if K is a boundary node. ! ! 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 IER > 0. ! ! NEAR,NEXT,DIST = Garbage. ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = -1 if N < 3 on input. ! IER = -2 if the first three nodes are ! collinear. ! IER = L if nodes L and M coincide for some ! M > L. The data structure represents ! a triangulation of nodes 1 to M-1 in ! this case. ! ! Modules required by TRMESH: ADDNOD, BDYADD, COVSPH, ! INSERT, INTADD, JRAND, ! LEFT, LSTPTR, STORE, SWAP, ! SWPTST, TRFIND ! ! Intrinsic function called by TRMESH: ABS ! !*********************************************************** ! INTEGER I, I0, J, K, LP, LPL, NEXTI, NN LOGICAL LEFT REAL D, D1, D2, D3 ! ! 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: ! K > 3 ! 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 ! NN = N IF (NN .LT. 3) THEN IER = -1 RETURN ENDIF ! ! 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 ! ELSEIF ( .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 ! ELSE ! ! The first three nodes are collinear. ! IER = -2 RETURN ENDIF ! ! Initialize LNEW and test for N = 3. ! LNEW = 7 IF (NN .EQ. 3) THEN IER = 0 RETURN ENDIF ! ! 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 I > 0) is the second, ! NEXT(L) (if L > 0) 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 1 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 .LE. D2 .AND. D1 .LE. D3) THEN NEAR(K) = 1 DIST(K) = D1 NEXT(K) = NEAR(1) NEAR(1) = K ELSEIF (D2 .LE. D1 .AND. D2 .LE. 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 ENDIF 1 CONTINUE ! ! Add the remaining nodes ! DO 6 K = 4,NN CALL ADDNOD (NEAR(K),K,X,Y,Z, LIST,LPTR,LEND, LNEW, IER) IF (IER .NE. 0) RETURN ! ! Remove K from the set of unprocessed nodes associated ! with NEAR(K). ! I = NEAR(K) IF (NEAR(I) .EQ. K) THEN NEAR(I) = NEXT(K) ELSE I = NEAR(I) 2 I0 = I I = NEXT(I0) IF (I .NE. K) GO TO 2 NEXT(I0) = NEXT(K) ENDIF NEAR(K) = 0 ! ! Loop on neighbors J of node K. ! LPL = LEND(K) LP = LPL 3 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) 4 IF (I .EQ. 0) GO TO 5 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 .LT. 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 .EQ. NEAR(J)) THEN NEAR(J) = NEXTI ELSE NEXT(I0) = NEXTI ENDIF NEXT(I) = NEAR(K) NEAR(K) = I ELSE I0 = I ENDIF ! ! Bottom of loop on I. ! I = NEXTI GO TO 4 ! ! Bottom of loop on neighbors J. ! 5 IF (LP .NE. LPL) GO TO 3 6 CONTINUE RETURN END SUBROUTINE TRPLOT (LUN,PLTSIZ,ELAT,ELON,A,N,X,Y,Z, & & LIST,LPTR,LEND,TITLE,NUMBR, IER) CHARACTER*(*) TITLE INTEGER LUN, N, LIST(*), LPTR(*), LEND(N), IER LOGICAL NUMBR REAL PLTSIZ, ELAT, ELON, A, X(N), Y(N), Z(N) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/16/98 ! ! This subroutine creates a level-2 Encapsulated Post- ! script (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. ! ! ! On input: ! ! LUN = 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. ! ! PLTSIZ = Plot size in inches. A circular window in ! the projection plane is mapped to a circu- ! lar 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 .LE. PLTSIZ .LE. 8.5. ! ! ELAT,ELON = 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 con- ! tains 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, respec- ! tively. ! ! A = Angular distance in degrees from E to the boun- ! dary 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 .LT. A ! .LE. 90. ! ! N = Number of nodes in the triangulation. N .GE. 3. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes (unit vectors). ! ! LIST,LPTR,LEND = Data structure defining the trian- ! gulation. Refer to Subroutine ! TRMESH. ! ! TITLE = Type CHARACTER variable or constant contain- ! ing 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 char- ! acters including the parentheses. ! ! NUMBR = Option indicator: If NUMBR = TRUE, the ! nodal indexes are plotted next to the nodes. ! ! Input parameters are not altered by this routine. ! ! On output: ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if LUN, PLTSIZ, or N is outside its ! valid range. ! IER = 2 if ELAT, ELON, or A is outside its ! valid range. ! IER = 3 if an error was encountered in writing ! to unit LUN. ! ! The values in the data statement below may be altered ! in order to modify various plotting options. ! ! Modules required by TRPLOT: None ! ! Intrinsic functions called by TRPLOT: ABS, ATAN, COS, ! NINT, REAL, SIN, ! SQRT ! !*********************************************************** ! INTEGER IPX1, IPX2, IPY1, IPY2, IR, LP, LPL, N0, N1 LOGICAL ANNOT REAL CF, CT, EX, EY, EZ, FSIZN, FSIZT, R11, R12, & & R21, R22, R23, SF, T, TX, TY, WR, WRS, X0, X1, & & Y0, Y1, Z0, Z1 ! DATA ANNOT/.TRUE./, FSIZN/10.0/, FSIZT/16.0/ ! ! 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) ! ! ! Test for invalid parameters. ! IF (LUN .LT. 0 .OR. LUN .GT. 99 .OR. & & PLTSIZ .LT. 1.0 .OR. PLTSIZ .GT. 8.5 .OR. & & N .LT. 3) & & GO TO 11 IF (ABS(ELAT) .GT. 90.0 .OR. ABS(ELON) .GT. 180.0 & & .OR. A .GT. 90.0) GO TO 12 ! ! Compute a conversion factor CF for degrees to radians ! and compute the window radius WR. ! CF = ATAN(1.0)/45.0 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.0*PLTSIZ) IPX1 = 306 - IR IPX2 = 306 + IR IPY1 = 396 - IR IPY2 = 396 + IR ! ! Output header comments. ! WRITE (LUN,100,ERR=13) IPX1, IPY1, IPX2, IPY2 100 FORMAT ('%!PS-Adobe-3.0 EPSF-3.0'/ & & '%%BoundingBox:',4I4/ & & '%%Title: Triangulation'/ & & '%%Creator: STRIPACK'/ & & '%%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.88*REAL(IR)) 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.0 WRITE (LUN,110,ERR=13) T WRITE (LUN,120,ERR=13) IR WRITE (LUN,130,ERR=13) 110 FORMAT (F12.6,' setlinewidth') 120 FORMAT ('306 396 ',I3,' 0 360 arc') 130 FORMAT ('stroke') ! ! Set up an affine mapping from the window box [-WR,WR] X ! [-WR,WR] to the viewport box. ! SF = REAL(IR)/WR TX = IPX1 + SF*WR TY = IPY1 + SF*WR WRITE (LUN,140,ERR=13) TX, TY, SF, SF 140 FORMAT (2F12.6,' translate'/ & & 2F12.6,' 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.0/SF WRITE (LUN,110,ERR=13) T ! ! Save the current graphics state, and set the clip path to ! the boundary of the window. ! WRITE (LUN,150,ERR=13) WRITE (LUN,160,ERR=13) WR WRITE (LUN,170,ERR=13) 150 FORMAT ('gsave') 160 FORMAT ('0 0 ',F12.6,' 0 360 arc') 170 FORMAT ('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 .NE. 0.0) THEN R11 = -EY/CT R12 = EX/CT ELSE R11 = 0.0 R12 = 1.0 ENDIF R21 = -EZ*R12 R22 = EZ*R11 R23 = CT ! ! Loop on visible nodes N0 that project to points (X0,Y0) in ! the window. ! DO 3 N0 = 1,N Z0 = EX*X(N0) + EY*Y(N0) + EZ*Z(N0) IF (Z0 .LT. 0.) GO TO 3 X0 = R11*X(N0) + R12*Y(N0) Y0 = R21*X(N0) + R22*Y(N0) + R23*Z(N0) IF (X0*X0 + Y0*Y0 .GT. WRS) GO TO 3 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. ! 1 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) IF (Z1 .LT. 0.) THEN ! ! 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. ! X1 = Z0*X1 - Z1*X0 Y1 = Z0*Y1 - Z1*Y0 T = SQRT(X1*X1+Y1*Y1) X1 = X1/T Y1 = Y1/T ENDIF ! ! If node N1 is in the window and N1 < N0, bypass edge ! N0->N1 (since edge N1->N0 has already been drawn). ! IF ( Z1 .GE. 0.0 .AND. X1*X1 + Y1*Y1 .LE. WRS & & .AND. N1 .LT. N0 ) GO TO 2 ! ! Add the edge to the path. ! WRITE (LUN,180,ERR=13) X0, Y0, X1, Y1 180 FORMAT (2F12.6,' moveto',2F12.6,' lineto') ! ! Bottom of loops. ! 2 IF (LP .NE. LPL) GO TO 1 3 CONTINUE ! ! Paint the path and restore the saved graphics state (with ! no clip path). ! WRITE (LUN,130,ERR=13) WRITE (LUN,190,ERR=13) 190 FORMAT ('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,200,ERR=13) T 200 FORMAT ('/Helvetica findfont'/ & & F12.6,' scalefont setfont') ! ! Loop on visible nodes N0 that project to points (X0,Y0) in ! the window. ! DO 4 N0 = 1,N IF (EX*X(N0) + EY*Y(N0) + EZ*Z(N0) .LT. 0.) & & GO TO 4 X0 = R11*X(N0) + R12*Y(N0) Y0 = R21*X(N0) + R22*Y(N0) + R23*Z(N0) IF (X0*X0 + Y0*Y0 .GT. WRS) GO TO 4 ! ! 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,210,ERR=13) X0, Y0 WRITE (LUN,220,ERR=13) N0 210 FORMAT (2F12.6,' moveto') 220 FORMAT ('(',I3,') show') 4 CONTINUE ENDIF ! ! Convert FSIZT from points to world coordinates, and output ! the commands to select a font and scale it. ! T = FSIZT/SF WRITE (LUN,200,ERR=13) T ! ! Display TITLE centered above the plot: ! Y0 = WR + 3.0*T WRITE (LUN,230,ERR=13) TITLE, Y0 230 FORMAT (A80/' stringwidth pop 2 div neg ',F12.6, ' moveto') WRITE (LUN,240,ERR=13) TITLE 240 FORMAT (A80/' show') IF (ANNOT) THEN ! ! Display the window center and radius below the plot. ! X0 = -WR Y0 = -WR - 50.0/SF WRITE (LUN,210,ERR=13) X0, Y0 WRITE (LUN,250,ERR=13) ELAT, ELON Y0 = Y0 - 2.0*T WRITE (LUN,210,ERR=13) X0, Y0 WRITE (LUN,260,ERR=13) A 250 FORMAT ('(Window center: ELAT = ',F7.2, & & ', ELON = ',F8.2,') show') 260 FORMAT ('(Angular extent: A = ',F5.2,') show') ENDIF ! ! Paint the path and output the showpage command and ! end-of-file indicator. ! WRITE (LUN,270,ERR=13) 270 FORMAT ('stroke'/ & & 'showpage'/ & & '%%EOF') ! ! HP's interpreters require a one-byte End-of-PostScript-Job ! indicator (to eliminate a timeout error message): ! ASCII 4. ! WRITE (LUN,280,ERR=13) CHAR(4) 280 FORMAT (A1) ! ! No error encountered. ! IER = 0 RETURN ! ! Invalid input parameter LUN, PLTSIZ, or N. ! 11 IER = 1 RETURN ! ! Invalid input parameter ELAT, ELON, or A. ! 12 IER = 2 RETURN ! ! Error writing to unit LUN. ! 13 IER = 3 RETURN END SUBROUTINE TRPRNT (N,X,Y,Z,IFLAG,LIST,LPTR,LEND,LOUT) INTEGER N, IFLAG, LIST(*), LPTR(*), LEND(N), LOUT REAL X(N), Y(N), Z(N) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/25/98 ! ! This subroutine prints the triangulation adjacency lists ! created by Subroutine 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. ! ! ! On input: ! ! N = Number of nodes in the triangulation. N .GE. 3 ! and N .LE. 9999. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes if IFLAG = 0, or ! (X and Y only) arrays of length N containing ! longitude and latitude, respectively, if ! IFLAG > 0, or unused dummy parameters if ! IFLAG < 0. ! ! IFLAG = Nodal coordinate option indicator: ! IFLAG = 0 if X, Y, and Z (assumed to contain ! Cartesian coordinates) are to be ! printed (to 6 decimal places). ! IFLAG > 0 if only X and Y (assumed to con- ! tain longitude and latitude) are ! to be printed (to 6 decimal ! places). ! IFLAG < 0 if only the adjacency lists are to ! be printed. ! ! LIST,LPTR,LEND = Data structure defining the trian- ! gulation. Refer to Subroutine ! TRMESH. ! ! LOUT = Logical unit for output. If LOUT is not in ! the range 0 to 99, output is written to ! logical unit 6. ! ! Input parameters are not altered by this routine. ! ! On output: ! ! The adjacency lists and nodal coordinates (as specified ! by IFLAG) are written to unit LOUT. ! ! Modules required by TRPRNT: None ! !*********************************************************** ! INTEGER I, INC, K, LP, LPL, LUN, NA, NABOR(400), NB, & & ND, NL, NLMAX, NMAX, NODE, NN, NT DATA NMAX/9999/, NLMAX/58/ ! ! 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 ! LUN = Logical unit for output (copy of LOUT) ! 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 addi- ! tional 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 ! NN = N LUN = LOUT IF (LUN .LT. 0 .OR. LUN .GT. 99) LUN = 6 ! ! Print a heading and test the range of N. ! WRITE (LUN,100) NN IF (NN .LT. 3 .OR. NN .GT. NMAX) THEN ! ! N is outside its valid range. ! WRITE (LUN,110) RETURN ENDIF ! ! Initialize NL (the number of lines printed on the current ! page) and NB (the number of boundary nodes encountered). ! NL = 6 NB = 0 IF (IFLAG .LT. 0) THEN ! ! Print LIST only. K is the number of neighbors of NODE ! that have been stored in NABOR. ! WRITE (LUN,101) DO 2 NODE = 1,NN LPL = LEND(NODE) LP = LPL K = 0 ! 1 K = K + 1 LP = LPTR(LP) ND = LIST(LP) NABOR(K) = ND IF (LP .NE. LPL) GO TO 1 IF (ND .LE. 0) THEN ! ! NODE is a boundary node. Correct the sign of the last ! neighbor, add 0 to the end of the list, and increment ! NB. ! NABOR(K) = -ND K = K + 1 NABOR(K) = 0 NB = NB + 1 ENDIF ! ! Increment NL and print the list of neighbors. ! INC = (K-1)/14 + 2 NL = NL + INC IF (NL .GT. NLMAX) THEN WRITE (LUN,108) NL = INC ENDIF WRITE (LUN,104) NODE, (NABOR(I), I = 1,K) IF (K .NE. 14) WRITE (LUN,107) 2 CONTINUE ELSEIF (IFLAG .GT. 0) THEN ! ! Print X (longitude), Y (latitude), and LIST. ! WRITE (LUN,102) DO 4 NODE = 1,NN LPL = LEND(NODE) LP = LPL K = 0 ! 3 K = K + 1 LP = LPTR(LP) ND = LIST(LP) NABOR(K) = ND IF (LP .NE. LPL) GO TO 3 IF (ND .LE. 0) THEN ! ! NODE is a boundary node. ! NABOR(K) = -ND K = K + 1 NABOR(K) = 0 NB = NB + 1 ENDIF ! ! Increment NL and print X, Y, and NABOR. ! INC = (K-1)/8 + 2 NL = NL + INC IF (NL .GT. NLMAX) THEN WRITE (LUN,108) NL = INC ENDIF WRITE (LUN,105) NODE, X(NODE), Y(NODE), (NABOR(I), I = 1,K) IF (K .NE. 8) WRITE (LUN,107) 4 CONTINUE ELSE ! ! Print X, Y, Z, and LIST. ! WRITE (LUN,103) DO 6 NODE = 1,NN LPL = LEND(NODE) LP = LPL K = 0 ! 5 K = K + 1 LP = LPTR(LP) ND = LIST(LP) NABOR(K) = ND IF (LP .NE. LPL) GO TO 5 IF (ND .LE. 0) THEN ! ! NODE is a boundary node. ! NABOR(K) = -ND K = K + 1 NABOR(K) = 0 NB = NB + 1 ENDIF ! ! Increment NL and print X, Y, Z, and NABOR. ! INC = (K-1)/5 + 2 NL = NL + INC IF (NL .GT. NLMAX) THEN WRITE (LUN,108) NL = INC ENDIF WRITE (LUN,106) NODE, X(NODE), Y(NODE), & & Z(NODE), (NABOR(I), I = 1,K) IF (K .NE. 5) WRITE (LUN,107) 6 CONTINUE ENDIF ! ! Print NB, NA, and NT (boundary nodes, arcs, and ! triangles). ! IF (NB .NE. 0) THEN NA = 3*NN - NB - 3 NT = 2*NN - NB - 2 ELSE NA = 3*NN - 6 NT = 2*NN - 4 ENDIF WRITE (LUN,109) NB, NA, NT RETURN ! ! Print formats: ! 100 FORMAT (///15X,'STRIPACK Triangulation Data ', & & 'Structure, N = ',I5//) 101 FORMAT (1X,'Node',31X,'Neighbors of Node'//) 102 FORMAT (1X,'Node',5X,'Longitude',6X,'Latitude', & & 18X,'Neighbors of Node'//) 103 FORMAT (1X,'Node',5X,'X(Node)',8X,'Y(Node)',8X, & & 'Z(Node)',11X,'Neighbors of Node'//) 104 FORMAT (1X,I4,4X,14I5/(1X,8X,14I5)) 105 FORMAT (1X,I4,2E15.6,4X,8I5/(1X,38X,8I5)) 106 FORMAT (1X,I4,3E15.6,4X,5I5/(1X,53X,5I5)) 107 FORMAT (1X) 108 FORMAT (///) 109 FORMAT (/1X,'NB = ',I4,' Boundary Nodes',5X, & & 'NA = ',I5,' Arcs',5X,'NT = ',I5, & & ' Triangles') 110 FORMAT (1X,10X,'*** N is outside its valid', & & ' range ***') END SUBROUTINE VRPLOT (LUN,PLTSIZ,ELAT,ELON,A,N,X,Y,Z, & & NT,LISTC,LPTR,LEND,XC,YC,ZC,TITLE, & & NUMBR, IER) CHARACTER*(*) TITLE INTEGER LUN, N, NT, LISTC(*), LPTR(*), LEND(N), IER LOGICAL NUMBR REAL PLTSIZ, ELAT, ELON, A, X(N), Y(N), Z(N), & & XC(NT), YC(NT), ZC(NT) ! !*********************************************************** ! ! From STRIPACK ! Robert J. Renka ! Dept. of Computer Science ! Univ. of North Texas ! renka@cs.unt.edu ! 07/16/98 ! ! This subroutine creates a level-2 Encapsulated Post- ! script (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 com- ! puted by Subroutine CRLIST. ! ! ! On input: ! ! LUN = 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. ! ! PLTSIZ = Plot size in inches. A circular window in ! the projection plane is mapped to a circu- ! lar 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 .LE. PLTSIZ .LE. 8.5. ! ! ELAT,ELON = 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 con- ! tains 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, respec- ! tively. ! ! A = Angular distance in degrees from E to the boun- ! dary of a circular window against which the ! Voronoi diagram is clipped. The projected win- ! dow 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 .LT. A .LE. 90. ! ! N = Number of nodes (Voronoi centers) and Voronoi ! regions. N .GE. 3. ! ! X,Y,Z = Arrays of length N containing the Cartesian ! coordinates of the nodes (unit vectors). ! ! NT = 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. ! ! LISTC = Array of length 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). ! ! LPTR = Array of length 3*NT = 6*N-12 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). ! ! LEND = Array of length N containing 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. ! ! XC,YC,ZC = Arrays of length NT containing the ! Cartesian coordinates of the triangle ! circumcenters (Voronoi vertices). ! XC(I)**2 + YC(I)**2 + ZC(I)**2 = 1. ! ! TITLE = Type CHARACTER variable or constant contain- ! ing 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 char- ! acters including the parentheses. ! ! NUMBR = Option indicator: If NUMBR = TRUE, the ! nodal indexes are plotted at the Voronoi ! region centers. ! ! Input parameters are not altered by this routine. ! ! On output: ! ! IER = Error indicator: ! IER = 0 if no errors were encountered. ! IER = 1 if LUN, PLTSIZ, N, or NT is outside ! its valid range. ! IER = 2 if ELAT, ELON, or A is outside its ! valid range. ! IER = 3 if an error was encountered in writing ! to unit LUN. ! ! Modules required by VRPLOT: None ! ! Intrinsic functions called by VRPLOT: ABS, ATAN, COS, ! NINT, REAL, SIN, ! SQRT ! !*********************************************************** ! INTEGER IPX1, IPX2, IPY1, IPY2, IR, KV1, KV2, LP, LPL, N0 LOGICAL ANNOT, IN1, IN2 REAL CF, CT, EX, EY, EZ, FSIZN, FSIZT, R11, R12, & & R21, R22, R23, SF, T, TX, TY, WR, WRS, X0, X1, & & X2, Y0, Y1, Y2, Z1, Z2 ! DATA ANNOT/.TRUE./, FSIZN/10.0/, FSIZT/16.0/ ! ! 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) ! ! ! Test for invalid parameters. ! IF (LUN .LT. 0 .OR. LUN .GT. 99 .OR. & & PLTSIZ .LT. 1.0 .OR. PLTSIZ .GT. 8.5 .OR. & & N .LT. 3 .OR. NT .NE. 2*N-4) & & GO TO 11 IF (ABS(ELAT) .GT. 90.0 .OR. ABS(ELON) .GT. 180.0 & & .OR. A .GT. 90.0) GO TO 12 ! ! Compute a conversion factor CF for degrees to radians ! and compute the window radius WR. ! CF = ATAN(1.0)/45.0 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.0*PLTSIZ) IPX1 = 306 - IR IPX2 = 306 + IR IPY1 = 396 - IR IPY2 = 396 + IR ! ! Output header comments. ! WRITE (LUN,100,ERR=13) IPX1, IPY1, IPX2, IPY2 100 FORMAT ('%!PS-Adobe-3.0 EPSF-3.0'/ & & '%%BoundingBox:',4I4/ & & '%%Title: Voronoi diagram'/ & & '%%Creator: STRIPACK'/ & & '%%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.88*REAL(IR)) 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.0 WRITE (LUN,110,ERR=13) T WRITE (LUN,120,ERR=13) IR WRITE (LUN,130,ERR=13) 110 FORMAT (F12.6,' setlinewidth') 120 FORMAT ('306 396 ',I3,' 0 360 arc') 130 FORMAT ('stroke') ! ! Set up an affine mapping from the window box [-WR,WR] X ! [-WR,WR] to the viewport box. ! SF = REAL(IR)/WR TX = IPX1 + SF*WR TY = IPY1 + SF*WR WRITE (LUN,140,ERR=13) TX, TY, SF, SF 140 FORMAT (2F12.6,' translate'/ & & 2F12.6,' 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.0/SF WRITE (LUN,110,ERR=13) T ! ! Save the current graphics state, and set the clip path to ! the boundary of the window. ! WRITE (LUN,150,ERR=13) WRITE (LUN,160,ERR=13) WR WRITE (LUN,170,ERR=13) 150 FORMAT ('gsave') 160 FORMAT ('0 0 ',F12.6,' 0 360 arc') 170 FORMAT ('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 .NE. 0.0) THEN R11 = -EY/CT R12 = EX/CT ELSE R11 = 0.0 R12 = 1.0 ENDIF R21 = -EZ*R12 R22 = EZ*R11 R23 = CT ! ! Loop on nodes (Voronoi centers) N0. ! LPL indexes the last neighbor of N0. ! DO 3 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 = Z2 .GE. 0. .AND. X2*X2 + Y2*Y2 .LE. WRS ! ! Loop on neighbors N1 of N0. For each triangulation edge ! N0-N1, KV1-KV2 is the corresponding Voronoi edge. ! LP = LPL 1 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 = Z2 .GE. 0. .AND. X2*X2 + Y2*Y2 .LE. WRS ! ! Add edge KV1-KV2 to the path iff both endpoints are inside ! the window and KV2 > KV1, or KV1 is inside and KV2 is ! outside (so that the edge is drawn only once). ! IF (.NOT. IN1 .OR. (IN2 .AND. KV2 .LE. KV1)) & & GO TO 2 IF (Z2 .LT. 0.) THEN ! ! 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. ! X2 = Z1*X2 - Z2*X1 Y2 = Z1*Y2 - Z2*Y1 T = SQRT(X2*X2+Y2*Y2) X2 = X2/T Y2 = Y2/T ENDIF WRITE (LUN,180,ERR=13) X1, Y1, X2, Y2 180 FORMAT (2F12.6,' moveto',2F12.6,' lineto') ! ! Bottom of loops. ! 2 IF (LP .NE. LPL) GO TO 1 3 CONTINUE ! ! Paint the path and restore the saved graphics state (with ! no clip path). ! WRITE (LUN,130,ERR=13) WRITE (LUN,190,ERR=13) 190 FORMAT ('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,200,ERR=13) T 200 FORMAT ('/Helvetica findfont'/ & & F12.6,' scalefont setfont') ! ! Loop on visible nodes N0 that project to points (X0,Y0) in ! the window. ! DO 4 N0 = 1,N IF (EX*X(N0) + EY*Y(N0) + EZ*Z(N0) .LT. 0.) & & GO TO 4 X0 = R11*X(N0) + R12*Y(N0) Y0 = R21*X(N0) + R22*Y(N0) + R23*Z(N0) IF (X0*X0 + Y0*Y0 .GT. WRS) GO TO 4 ! ! Move to (X0,Y0), and draw the label N0 with the origin ! of the first character at (X0,Y0). ! WRITE (LUN,210,ERR=13) X0, Y0 WRITE (LUN,220,ERR=13) N0 210 FORMAT (2F12.6,' moveto') 220 FORMAT ('(',I3,') show') 4 CONTINUE ENDIF ! ! Convert FSIZT from points to world coordinates, and output ! the commands to select a font and scale it. ! T = FSIZT/SF WRITE (LUN,200,ERR=13) T ! ! Display TITLE centered above the plot: ! Y0 = WR + 3.0*T WRITE (LUN,230,ERR=13) TITLE, Y0 230 FORMAT (A80/' stringwidth pop 2 div neg ',F12.6, & & ' moveto') WRITE (LUN,240,ERR=13) TITLE 240 FORMAT (A80/' show') IF (ANNOT) THEN ! ! Display the window center and radius below the plot. ! X0 = -WR Y0 = -WR - 50.0/SF WRITE (LUN,210,ERR=13) X0, Y0 WRITE (LUN,250,ERR=13) ELAT, ELON Y0 = Y0 - 2.0*T WRITE (LUN,210,ERR=13) X0, Y0 WRITE (LUN,260,ERR=13) A 250 FORMAT ('(Window center: ELAT = ',F7.2, & & ', ELON = ',F8.2,') show') 260 FORMAT ('(Angular extent: A = ',F5.2,') show') ENDIF ! ! Paint the path and output the showpage command and ! end-of-file indicator. ! WRITE (LUN,270,ERR=13) 270 FORMAT ('stroke'/ & & 'showpage'/ & & '%%EOF') ! ! HP's interpreters require a one-byte End-of-PostScript-Job ! indicator (to eliminate a timeout error message): ! ASCII 4. ! WRITE (LUN,280,ERR=13) CHAR(4) 280 FORMAT (A1) ! ! No error encountered. ! IER = 0 RETURN ! ! Invalid input parameter LUN, PLTSIZ, N, or NT. ! 11 IER = 1 RETURN ! ! Invalid input parameter ELAT, ELON, or A. ! 12 IER = 2 RETURN ! ! Error writing to unit LUN. ! 13 IER = 3 RETURN END