= det(N0,N1,N2) .GE. 0.
C
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(*)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/15/96
C
C This function returns the index (LIST pointer) of NB in
C the adjacency list for N0, where LPL = LEND(N0).
C
C This function is identical to the similarly named
C function in TRIPACK.
C
C
C On input:
C
C LPL = LEND(N0)
C
C NB = Index of the node whose pointer is to be re-
C turned. NB must be connected to N0.
C
C LIST,LPTR = Data structure defining the triangula-
C tion. Refer to Subroutine TRMESH.
C
C Input parameters are not altered by this function.
C
C On output:
C
C LSTPTR = Pointer such that LIST(LSTPTR) = NB or
C LIST(LSTPTR) = -NB, unless NB is not a
C neighbor of N0, in which case LSTPTR = LPL.
C
C Modules required by LSTPTR: None
C
C***********************************************************
C
INTEGER LP, ND
C
C Local parameters:
C
C LP = LIST pointer
C ND = Nodal index
C
LP = LPTR(LPL)
1 ND = LIST(LP)
IF (ND .EQ. NB) GO TO 2
LP = LPTR(LP)
IF (LP .NE. LPL) GO TO 1
C
2 LSTPTR = LP
RETURN
END
INTEGER FUNCTION NBCNT (LPL,LPTR)
INTEGER LPL, LPTR(*)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/15/96
C
C This function returns the number of neighbors of a node
C N0 in a triangulation created by Subroutine TRMESH.
C
C This function is identical to the similarly named
C function in TRIPACK.
C
C
C On input:
C
C LPL = LIST pointer to the last neighbor of N0 --
C LPL = LEND(N0).
C
C LPTR = Array of pointers associated with LIST.
C
C Input parameters are not altered by this function.
C
C On output:
C
C NBCNT = Number of neighbors of N0.
C
C Modules required by NBCNT: None
C
C***********************************************************
C
INTEGER K, LP
C
C Local parameters:
C
C K = Counter for computing the number of neighbors
C LP = LIST pointer
C
LP = LPL
K = 1
C
1 LP = LPTR(LP)
IF (LP .EQ. LPL) GO TO 2
K = K + 1
GO TO 1
C
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
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/28/98
C
C Given a point P on the surface of the unit sphere and a
C Delaunay triangulation created by Subroutine TRMESH, this
C function returns the index of the nearest triangulation
C node to P.
C
C The algorithm consists of implicitly adding P to the
C triangulation, finding the nearest neighbor to P, and
C implicitly deleting P from the triangulation. Thus, it
C is based on the fact that, if P is a node in a Delaunay
C triangulation, the nearest node to P is a neighbor of P.
C
C
C On input:
C
C P = Array of length 3 containing the Cartesian coor-
C dinates of the point P to be located relative to
C the triangulation. It is assumed without a test
C that P(1)**2 + P(2)**2 + P(3)**2 = 1.
C
C IST = Index of a node at which TRFIND begins the
C search. Search time depends on the proximity
C of this node to P.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y,Z = Arrays of length N containing the Cartesian
C coordinates of the nodes.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to TRMESH.
C
C Input parameters are not altered by this function.
C
C On output:
C
C NEARND = Nodal index of the nearest node to P, or 0
C if N < 3 or the triangulation data struc-
C ture is invalid.
C
C AL = Arc length (angular distance in radians) be-
C tween P and NEARND unless NEARND = 0.
C
C Note that the number of candidates for NEARND
C (neighbors of P) is limited to LMAX defined in
C the PARAMETER statement below.
C
C Modules required by NEARND: JRAND, LSTPTR, TRFIND, STORE
C
C Intrinsic functions called by NEARND: ABS, ACOS
C
C***********************************************************
C
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
C
C Local parameters:
C
C B1,B2,B3 = Unnormalized barycentric coordinates returned
C by TRFIND
C DS1 = (Negative cosine of the) distance from P to N1
C DSR = (Negative cosine of the) distance from P to NR
C DX1,..DZ3 = Components of vectors used by the swap test
C I1,I2,I3 = Nodal indexes of a triangle containing P, or
C the rightmost (I1) and leftmost (I2) visible
C boundary nodes as viewed from P
C L = Length of LISTP/LPTRP and number of neighbors
C of P
C LMAX = Maximum value of L
C LISTP = Indexes of the neighbors of P
C LPTRP = Array of pointers in 1-1 correspondence with
C LISTP elements
C LP = LIST pointer to a neighbor of N1 and LISTP
C pointer
C LP1,LP2 = LISTP indexes (pointers)
C LPL = Pointer to the last neighbor of N1
C N1 = Index of a node visible from P
C N2 = Index of an endpoint of an arc opposite P
C N3 = Index of the node opposite N1->N2
C NN = Local copy of N
C NR = Index of a candidate for the nearest node to P
C NST = Index of the node at which TRFIND begins the
C search
C
C
C Store local parameters and test for N invalid.
C
NN = N
IF (NN .LT. 3) GO TO 6
NST = IST
IF (NST .LT. 1 .OR. NST .GT. NN) NST = 1
C
C Find a triangle (I1,I2,I3) containing P, or the rightmost
C (I1) and leftmost (I2) visible boundary nodes as viewed
C from P.
C
CALL TRFIND (NST,P,N,X,Y,Z,LIST,LPTR,LEND, B1,B2,B3,
. I1,I2,I3)
C
C Test for collinear nodes.
C
IF (I1 .EQ. 0) GO TO 6
C
C Store the linked list of 'neighbors' of P in LISTP and
C LPTRP. I1 is the first neighbor, and 0 is stored as
C the last neighbor if P is not contained in a triangle.
C L is the length of LISTP and LPTRP, and is limited to
C LMAX.
C
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
C
C Loop on the ordered sequence of visible boundary nodes
C N1 from I1 to I2.
C
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
C
C Initialize variables for a loop on arcs N1-N2 opposite P
C in which new 'neighbors' are 'swapped' in. N1 follows
C N2 as a neighbor of P, and LP1 and LP2 are the LISTP
C indexes of N1 and N2.
C
LP2 = 1
N2 = I1
LP1 = LPTRP(1)
N1 = LISTP(LP1)
C
C Begin loop: find the node N3 opposite N1->N2.
C
2 LP = LSTPTR(LEND(N1),N2,LIST,LPTR)
IF (LIST(LP) .LT. 0) GO TO 3
LP = LPTR(LP)
N3 = ABS(LIST(LP))
C
C Swap test: Exit the loop if L = LMAX.
C
IF (L .EQ. LMAX) GO TO 4
DX1 = X(N1) - P(1)
DY1 = Y(N1) - P(2)
DZ1 = Z(N1) - P(3)
C
DX2 = X(N2) - P(1)
DY2 = Y(N2) - P(2)
DZ2 = Z(N2) - P(3)
C
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
C
C Swap: Insert N3 following N2 in the adjacency list for P.
C The two new arcs opposite P must be tested.
C
L = L+1
LPTRP(LP2) = L
LISTP(L) = N3
LPTRP(L) = LP1
LP1 = L
N1 = N3
GO TO 2
C
C No swap: Advance to the next arc and test for termination
C on N1 = I1 (LP1 = 1) or N1 followed by 0.
C
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
C
C Set NR and DSR to the index of the nearest node to P and
C an increasing function (negative cosine) of its distance
C from P, respectively.
C
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
C
C Invalid input.
C
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(*)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/30/98
C
C Given a set of NA triangulation arcs, this subroutine
C optimizes the portion of the triangulation consisting of
C the quadrilaterals (pairs of adjacent triangles) which
C have the arcs as diagonals by applying the circumcircle
C test and appropriate swaps to the arcs.
C
C An iteration consists of applying the swap test and
C swaps to all NA arcs in the order in which they are
C stored. The iteration is repeated until no swap occurs
C or NIT iterations have been performed. The bound on the
C number of iterations may be necessary to prevent an
C infinite loop caused by cycling (reversing the effect of a
C previous swap) due to floating point inaccuracy when four
C or more nodes are nearly cocircular.
C
C
C On input:
C
C X,Y,Z = Arrays containing the nodal coordinates.
C
C NA = Number of arcs in the set. NA .GE. 0.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C NIT = Maximum number of iterations to be performed.
C NIT = 4*NA should be sufficient. NIT .GE. 1.
C
C IWK = Integer array dimensioned 2 by NA containing
C the nodal indexes of the arc endpoints (pairs
C of endpoints are stored in columns).
C
C On output:
C
C LIST,LPTR,LEND = Updated triangulation data struc-
C ture reflecting the swaps.
C
C NIT = Number of iterations performed.
C
C IWK = Endpoint indexes of the new set of arcs
C reflecting the swaps.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if a swap occurred on the last of
C MAXIT iterations, where MAXIT is the
C value of NIT on input. The new set
C of arcs is not necessarily optimal
C in this case.
C IER = 2 if NA < 0 or NIT < 1 on input.
C IER = 3 if IWK(2,I) is not a neighbor of
C IWK(1,I) for some I in the range 1
C to NA. A swap may have occurred in
C this case.
C IER = 4 if a zero pointer was returned by
C Subroutine SWAP.
C
C Modules required by OPTIM: LSTPTR, SWAP, SWPTST
C
C Intrinsic function called by OPTIM: ABS
C
C***********************************************************
C
INTEGER I, IO1, IO2, ITER, LP, LP21, LPL, LPP, MAXIT,
. N1, N2, NNA
LOGICAL SWPTST
LOGICAL SWP
C
C Local parameters:
C
C I = Column index for IWK
C IO1,IO2 = Nodal indexes of the endpoints of an arc in IWK
C ITER = Iteration count
C LP = LIST pointer
C LP21 = Parameter returned by SWAP (not used)
C LPL = Pointer to the last neighbor of IO1
C LPP = Pointer to the node preceding IO2 as a neighbor
C of IO1
C MAXIT = Input value of NIT
C N1,N2 = Nodes opposite IO1->IO2 and IO2->IO1,
C respectively
C NNA = Local copy of NA
C SWP = Flag set to TRUE iff a swap occurs in the
C optimization loop
C
NNA = NA
MAXIT = NIT
IF (NNA .LT. 0 .OR. MAXIT .LT. 1) GO TO 7
C
C Initialize iteration count ITER and test for NA = 0.
C
ITER = 0
IF (NNA .EQ. 0) GO TO 5
C
C Top of loop --
C SWP = TRUE iff a swap occurred in the current iteration.
C
1 IF (ITER .EQ. MAXIT) GO TO 6
ITER = ITER + 1
SWP = .FALSE.
C
C Inner loop on arcs IO1-IO2 --
C
DO 4 I = 1,NNA
IO1 = IWK(1,I)
IO2 = IWK(2,I)
C
C Set N1 and N2 to the nodes opposite IO1->IO2 and
C IO2->IO1, respectively. Determine the following:
C
C LPL = pointer to the last neighbor of IO1,
C LP = pointer to IO2 as a neighbor of IO1, and
C LPP = pointer to the node N2 preceding IO2.
C
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
C
C IO2 should be the last neighbor of IO1. Test for no
C arc and bypass the swap test if IO1 is a boundary
C node.
C
IF (ABS(LIST(LP)) .NE. IO2) GO TO 8
IF (LIST(LP) .LT. 0) GO TO 4
C
C Store N1 and N2, or bypass the swap test if IO1 is a
C boundary node and IO2 is its first neighbor.
C
3 N2 = LIST(LPP)
IF (N2 .LT. 0) GO TO 4
LP = LPTR(LP)
N1 = ABS(LIST(LP))
C
C Test IO1-IO2 for a swap, and update IWK if necessary.
C
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
C
C Successful termination.
C
5 NIT = ITER
IER = 0
RETURN
C
C MAXIT iterations performed without convergence.
C
6 NIT = MAXIT
IER = 1
RETURN
C
C Invalid input parameter.
C
7 NIT = 0
IER = 2
RETURN
C
C IO2 is not a neighbor of IO1.
C
8 NIT = ITER
IER = 3
RETURN
C
C Zero pointer returned by SWAP.
C
9 NIT = ITER
IER = 4
RETURN
END
SUBROUTINE SCOORD (PX,PY,PZ, PLAT,PLON,PNRM)
REAL PX, PY, PZ, PLAT, PLON, PNRM
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 08/27/90
C
C This subroutine converts a point P from Cartesian coor-
C dinates to spherical coordinates.
C
C
C On input:
C
C PX,PY,PZ = Cartesian coordinates of P.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C PLAT = Latitude of P in the range -PI/2 to PI/2, or
C 0 if PNRM = 0. PLAT should be scaled by
C 180/PI to obtain the value in degrees.
C
C PLON = Longitude of P in the range -PI to PI, or 0
C if P lies on the Z-axis. PLON should be
C scaled by 180/PI to obtain the value in
C degrees.
C
C PNRM = Magnitude (Euclidean norm) of P.
C
C Modules required by SCOORD: None
C
C Intrinsic functions called by SCOORD: ASIN, ATAN2, SQRT
C
C***********************************************************
C
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
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 05/09/92
C
C This function forces its argument X to be stored in a
C memory location, thus providing a means of determining
C floating point number characteristics (such as the machine
C precision) when it is necessary to avoid computation in
C high precision registers.
C
C
C On input:
C
C X = Value to be stored.
C
C X is not altered by this function.
C
C On output:
C
C STORE = Value of X after it has been stored and
C possibly truncated or rounded to the single
C precision word length.
C
C Modules required by STORE: None
C
C***********************************************************
C
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
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 06/22/98
C
C Given a triangulation of a set of points on the unit
C sphere, this subroutine replaces a diagonal arc in a
C strictly convex quadrilateral (defined by a pair of adja-
C cent triangles) with the other diagonal. Equivalently, a
C pair of adjacent triangles is replaced by another pair
C having the same union.
C
C
C On input:
C
C IN1,IN2,IO1,IO2 = Nodal indexes of the vertices of
C the quadrilateral. IO1-IO2 is re-
C placed by IN1-IN2. (IO1,IO2,IN1)
C and (IO2,IO1,IN2) must be trian-
C gles on input.
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C On output:
C
C LIST,LPTR,LEND = Data structure updated with the
C swap -- triangles (IO1,IO2,IN1) and
C (IO2,IO1,IN2) are replaced by
C (IN1,IN2,IO2) and (IN2,IN1,IO1)
C unless LP21 = 0.
C
C LP21 = Index of IN1 as a neighbor of IN2 after the
C swap is performed unless IN1 and IN2 are
C adjacent on input, in which case LP21 = 0.
C
C Module required by SWAP: LSTPTR
C
C Intrinsic function called by SWAP: ABS
C
C***********************************************************
C
INTEGER LSTPTR
INTEGER LP, LPH, LPSAV
C
C Local parameters:
C
C LP,LPH,LPSAV = LIST pointers
C
C
C Test for IN1 and IN2 adjacent.
C
LP = LSTPTR(LEND(IN1),IN2,LIST,LPTR)
IF (ABS(LIST(LP)) .EQ. IN2) THEN
LP21 = 0
RETURN
ENDIF
C
C Delete IO2 as a neighbor of IO1.
C
LP = LSTPTR(LEND(IO1),IN2,LIST,LPTR)
LPH = LPTR(LP)
LPTR(LP) = LPTR(LPH)
C
C If IO2 is the last neighbor of IO1, make IN2 the
C last neighbor.
C
IF (LEND(IO1) .EQ. LPH) LEND(IO1) = LP
C
C Insert IN2 as a neighbor of IN1 following IO1
C using the hole created above.
C
LP = LSTPTR(LEND(IN1),IO1,LIST,LPTR)
LPSAV = LPTR(LP)
LPTR(LP) = LPH
LIST(LPH) = IN2
LPTR(LPH) = LPSAV
C
C Delete IO1 as a neighbor of IO2.
C
LP = LSTPTR(LEND(IO2),IN1,LIST,LPTR)
LPH = LPTR(LP)
LPTR(LP) = LPTR(LPH)
C
C If IO1 is the last neighbor of IO2, make IN1 the
C last neighbor.
C
IF (LEND(IO2) .EQ. LPH) LEND(IO2) = LP
C
C Insert IN1 as a neighbor of IN2 following IO2.
C
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(*)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 03/29/91
C
C This function decides whether or not to replace a
C diagonal arc in a quadrilateral with the other diagonal.
C The decision will be to swap (SWPTST = TRUE) if and only
C if N4 lies above the plane (in the half-space not contain-
C ing the origin) defined by (N1,N2,N3), or equivalently, if
C the projection of N4 onto this plane is interior to the
C circumcircle of (N1,N2,N3). The decision will be for no
C swap if the quadrilateral is not strictly convex.
C
C
C On input:
C
C N1,N2,N3,N4 = Indexes of the four nodes defining the
C quadrilateral with N1 adjacent to N2,
C and (N1,N2,N3) in counterclockwise
C order. The arc connecting N1 to N2
C should be replaced by an arc connec-
C ting N3 to N4 if SWPTST = TRUE. Refer
C to Subroutine SWAP.
C
C X,Y,Z = Arrays of length N containing the Cartesian
C coordinates of the nodes. (X(I),Y(I),Z(I))
C define node I for I = N1, N2, N3, and N4.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C SWPTST = TRUE if and only if the arc connecting N1
C and N2 should be swapped for an arc con-
C necting N3 and N4.
C
C Modules required by SWPTST: None
C
C***********************************************************
C
REAL DX1, DX2, DX3, DY1, DY2, DY3, DZ1, DZ2, DZ3,
. X4, Y4, Z4
C
C Local parameters:
C
C DX1,DY1,DZ1 = Coordinates of N4->N1
C DX2,DY2,DZ2 = Coordinates of N4->N2
C DX3,DY3,DZ3 = Coordinates of N4->N3
C X4,Y4,Z4 = Coordinates of N4
C
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
C
C N4 lies above the plane of (N1,N2,N3) iff N3 lies above
C the plane of (N2,N1,N4) iff Det(N3-N4,N2-N4,N1-N4) =
C (N3-N4,N2-N4 X N1-N4) > 0.
C
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)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 04/08/90
C
C This subroutine transforms spherical coordinates into
C Cartesian coordinates on the unit sphere for input to
C Subroutine TRMESH. Storage for X and Y may coincide with
C storage for RLAT and RLON if the latter need not be saved.
C
C
C On input:
C
C N = Number of nodes (points on the unit sphere)
C whose coordinates are to be transformed.
C
C RLAT = Array of length N containing latitudinal
C coordinates of the nodes in radians.
C
C RLON = Array of length N containing longitudinal
C coordinates of the nodes in radians.
C
C The above parameters are not altered by this routine.
C
C X,Y,Z = Arrays of length at least N.
C
C On output:
C
C X,Y,Z = Cartesian coordinates in the range -1 to 1.
C X(I)**2 + Y(I)**2 + Z(I)**2 = 1 for I = 1
C to N.
C
C Modules required by TRANS: None
C
C Intrinsic functions called by TRANS: COS, SIN
C
C***********************************************************
C
INTEGER I, NN
REAL COSPHI, PHI, THETA
C
C Local parameters:
C
C COSPHI = cos(PHI)
C I = DO-loop index
C NN = Local copy of N
C PHI = Latitude
C THETA = Longitude
C
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
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 11/30/99
C
C This subroutine locates a point P relative to a triangu-
C lation created by Subroutine TRMESH. If P is contained in
C a triangle, the three vertex indexes and barycentric coor-
C dinates are returned. Otherwise, the indexes of the
C visible boundary nodes are returned.
C
C
C On input:
C
C NST = Index of a node at which TRFIND begins its
C search. Search time depends on the proximity
C of this node to P.
C
C P = Array of length 3 containing the x, y, and z
C coordinates (in that order) of the point P to be
C located.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y,Z = Arrays of length N containing the Cartesian
C coordinates of the triangulation nodes (unit
C vectors). (X(I),Y(I),Z(I)) defines node I
C for I = 1 to N.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C B1,B2,B3 = Unnormalized barycentric coordinates of
C the central projection of P onto the un-
C derlying planar triangle if P is in the
C convex hull of the nodes. These parame-
C ters are not altered if I1 = 0.
C
C I1,I2,I3 = Counterclockwise-ordered vertex indexes
C of a triangle containing P if P is con-
C tained in a triangle. If P is not in the
C convex hull of the nodes, I1 and I2 are
C the rightmost and leftmost (boundary)
C nodes that are visible from P, and
C I3 = 0. (If all boundary nodes are vis-
C ible from P, then I1 and I2 coincide.)
C I1 = I2 = I3 = 0 if P and all of the
C nodes are coplanar (lie on a common great
C circle.
C
C Modules required by TRFIND: JRAND, LSTPTR, STORE
C
C Intrinsic function called by TRFIND: ABS
C
C***********************************************************
C
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
C
SAVE IX, IY, IZ
DATA IX/1/, IY/2/, IZ/3/
C
C Local parameters:
C
C EPS = Machine precision
C IX,IY,IZ = Integer seeds for JRAND
C LP = LIST pointer
C N0,N1,N2 = Nodes in counterclockwise order defining a
C cone (with vertex N0) containing P, or end-
C points of a boundary edge such that P Right
C N1->N2
C N1S,N2S = Initially-determined values of N1 and N2
C N3,N4 = Nodes opposite N1->N2 and N2->N1, respectively
C NEXT = Candidate for I1 or I2 when P is exterior
C NF,NL = First and last neighbors of N0, or first
C (rightmost) and last (leftmost) nodes
C visible from P when P is exterior to the
C triangulation
C PTN1 = Scalar product
C PTN2 = Scalar product

C Q = (N2 X N1) X N2 or N1 X (N2 X N1) -- used in
C the boundary traversal when P is exterior
C S12 = Scalar product
C TOL = Tolerance (multiple of EPS) defining an upper
C bound on the magnitude of a negative bary-
C centric coordinate (B1 or B2) for P in a
C triangle -- used to avoid an infinite number
C of restarts with 0 <= B3 < EPS and B1 < 0 or
C B2 < 0 but small in magnitude
C XP,YP,ZP = Local variables containing P(1), P(2), and P(3)
C X0,Y0,Z0 = Dummy arguments for DET
C X1,Y1,Z1 = Dummy arguments for DET
C X2,Y2,Z2 = Dummy arguments for DET
C
C Statement function:
C
C DET(X1,...,Z0) .GE. 0 if and only if (X0,Y0,Z0) is in the
C (closed) left hemisphere defined by
C the plane containing (0,0,0),
C (X1,Y1,Z1), and (X2,Y2,Z2), where
C left is defined relative to an ob-
C server at (X1,Y1,Z1) facing
C (X2,Y2,Z2).
C
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)
C
C Initialize variables.
C
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 )
C
C Compute the relative machine precision EPS and TOL.
C
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
C
C Set NF and NL to the first and last neighbors of N0, and
C initialize N1 = NF.
C
2 LP = LEND(N0)
NL = LIST(LP)
LP = LPTR(LP)
NF = LIST(LP)
N1 = NF
C
C Find a pair of adjacent neighbors N1,N2 of N0 that define
C a wedge containing P: P LEFT N0->N1 and P RIGHT N0->N2.
C
IF (NL .GT. 0) THEN
C
C N0 is an interior node. Find N1.
C
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
C
C N0 is a boundary node. Test for P exterior.
C
NL = -NL
IF ( DET(X(N0),Y(N0),Z(N0),X(NF),Y(NF),Z(NF),
. XP,YP,ZP) .LT. 0. ) THEN
C
C P is to the right of the boundary edge N0->NF.
C
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
C
C P is to the right of the boundary edge NL->N0.
C
N1 = NL
N2 = N0
GO TO 9
ENDIF
ENDIF
C
C P is to the left of arcs N0->N1 and NL->N0. Set N2 to the
C next neighbor of N0 (following N1).
C
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
C
C P is left of or on arcs N0->NB for all neighbors NB
C of N0. Test for P = +/-N0.
C
IF (STORE(ABS(X(N0)*XP + Y(N0)*YP + Z(N0)*ZP))
. .LT. 1.0-4.0*EPS) THEN
C
C All points are collinear iff P Left NB->N0 for all
C neighbors NB of N0. Search the neighbors of N0.
C Note: N1 = NL and LP points to NL.
C
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
C
C P is to the right of N1->N0, or P = +/-N0. Set N0 to N1
C and start over.
C
N0 = N1
GO TO 2
C
C P is between arcs N0->N1 and N0->NF.
C
6 N2 = NF
C
C P is contained in a wedge defined by geodesics N0-N1 and
C N0-N2, where N1 is adjacent to N2. Save N1 and N2 to
C test for cycling.
C
7 N3 = N0
N1S = N1
N2S = N2
C
C Top of edge-hopping loop:
C
8 B3 = DET(X(N1),Y(N1),Z(N1),X(N2),Y(N2),Z(N2),XP,YP,ZP)
IF (B3 .LT. 0.) THEN
C
C Set N4 to the first neighbor of N2 following N1 (the
C node opposite N2->N1) unless N1->N2 is a boundary arc.
C
LP = LSTPTR(LEND(N2),N1,LIST,LPTR)
IF (LIST(LP) .LT. 0) GO TO 9
LP = LPTR(LP)
N4 = ABS(LIST(LP))
C
C Define a new arc N1->N2 which intersects the geodesic
C N0-P.
C
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
C
C The starting node N0 or edge N1-N2 was encountered
C again, implying a cycle (infinite loop). Restart
C with N0 randomly selected.
C
N0 = JRAND(N, IX,IY,IZ )
GO TO 2
ENDIF
C
C P is in (N1,N2,N3) unless N0, N1, N2, and P are collinear
C or P is close to -N0.
C
IF (B3 .GE. EPS) THEN
C
C B3 .NE. 0.
C
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
C
C Restart with N0 randomly selected.
C
N0 = JRAND(N, IX,IY,IZ )
GO TO 2
ENDIF
ELSE
C
C B3 = 0 and thus P lies on N1->N2. Compute
C B1 = Det(P,N2 X N1,N2) and B2 = Det(P,N1,N2 X N1).
C
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
C
C Restart with N0 randomly selected.
C
N0 = JRAND(N, IX,IY,IZ )
GO TO 2
ENDIF
ENDIF
C
C P is in (N1,N2,N3).
C
I1 = N1
I2 = N2
I3 = N3
IF (B1 .LT. 0.0) B1 = 0.0
IF (B2 .LT. 0.0) B2 = 0.0
RETURN
C
C P Right N1->N2, where N1->N2 is a boundary edge.
C Save N1 and N2, and set NL = 0 to indicate that
C NL has not yet been found.
C
9 N1S = N1
N2S = N2
NL = 0
C
C Counterclockwise Boundary Traversal:
C
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
C
C N2 is the rightmost visible node if P Forward N2->N1
C or NEXT Forward N2->N1. Set Q to (N2 X N1) X N2.
C
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
C
C N1, N2, NEXT, and P are nearly collinear, and N2 is
C the leftmost visible node.
C
NL = N2
ENDIF
C
C Bottom of counterclockwise loop:
C
N1 = N2
N2 = NEXT
IF (N2 .NE. N1S) GO TO 10
C
C All boundary nodes are visible from P.
C
I1 = N1S
I2 = N1S
I3 = 0
RETURN
C
C N2 is the rightmost visible node.
C
11 NF = N2
IF (NL .EQ. 0) THEN
C
C Restore initial values of N1 and N2, and begin the search
C for the leftmost visible node.
C
N2 = N2S
N1 = N1S
C
C Clockwise Boundary Traversal:
C
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
C
C N1 is the leftmost visible node if P or NEXT is
C forward of N1->N2. Compute Q = N1 X (N2 X N1).
C
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
C
C P, NEXT, N1, and N2 are nearly collinear and N1 is the
C rightmost visible node.
C
NF = N1
ENDIF
C
C Bottom of clockwise loop:
C
N2 = N1
N1 = NEXT
IF (N1 .NE. N1S) GO TO 12
C
C All boundary nodes are visible from P.
C
I1 = N1
I2 = N1
I3 = 0
RETURN
C
C N1 is the leftmost visible node.
C
13 NL = N1
ENDIF
C
C NF and NL have been found.
C
I1 = NF
I2 = NL
I3 = 0
RETURN
C
C All points are collinear (coplanar).
C
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
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/20/96
C
C This subroutine converts a triangulation data structure
C from the linked list created by Subroutine TRMESH to a
C triangle list.
C
C On input:
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C LIST,LPTR,LEND = Linked list data structure defin-
C ing the triangulation. Refer to
C Subroutine TRMESH.
C
C NROW = Number of rows (entries per triangle) re-
C served for the triangle list LTRI. The value
C must be 6 if only the vertex indexes and
C neighboring triangle indexes are to be
C stored, or 9 if arc indexes are also to be
C assigned and stored. Refer to LTRI.
C
C The above parameters are not altered by this routine.
C
C LTRI = Integer array of length at least NROW*NT,
C where NT is at most 2N-4. (A sufficient
C length is 12N if NROW=6 or 18N if NROW=9.)
C
C On output:
C
C NT = Number of triangles in the triangulation unless
C IER .NE. 0, in which case NT = 0. NT = 2N-NB-2
C if NB .GE. 3 or 2N-4 if NB = 0, where NB is the
C number of boundary nodes.
C
C LTRI = NROW by NT array whose J-th column contains
C the vertex nodal indexes (first three rows),
C neighboring triangle indexes (second three
C rows), and, if NROW = 9, arc indexes (last
C three rows) associated with triangle J for
C J = 1,...,NT. The vertices are ordered
C counterclockwise with the first vertex taken
C to be the one with smallest index. Thus,
C LTRI(2,J) and LTRI(3,J) are larger than
C LTRI(1,J) and index adjacent neighbors of
C node LTRI(1,J). For I = 1,2,3, LTRI(I+3,J)
C and LTRI(I+6,J) index the triangle and arc,
C respectively, which are opposite (not shared
C by) node LTRI(I,J), with LTRI(I+3,J) = 0 if
C LTRI(I+6,J) indexes a boundary arc. Vertex
C indexes range from 1 to N, triangle indexes
C from 0 to NT, and, if included, arc indexes
C from 1 to NA, where NA = 3N-NB-3 if NB .GE. 3
C or 3N-6 if NB = 0. The triangles are or-
C dered on first (smallest) vertex indexes.
C
C IER = Error indicator.
C IER = 0 if no errors were encountered.
C IER = 1 if N or NROW is outside its valid
C range on input.
C IER = 2 if the triangulation data structure
C (LIST,LPTR,LEND) is invalid. Note,
C however, that these arrays are not
C completely tested for validity.
C
C Modules required by TRLIST: None
C
C Intrinsic function called by TRLIST: ABS
C
C***********************************************************
C
INTEGER I, I1, I2, I3, ISV, J, KA, KN, KT, LP, LP2,
. LPL, LPLN1, N1, N2, N3, NM2
LOGICAL ARCS
C
C Local parameters:
C
C ARCS = Logical variable with value TRUE iff are
C indexes are to be stored
C I,J = LTRI row indexes (1 to 3) associated with
C triangles KT and KN, respectively
C I1,I2,I3 = Nodal indexes of triangle KN
C ISV = Variable used to permute indexes I1,I2,I3
C KA = Arc index and number of currently stored arcs
C KN = Index of the triangle that shares arc I1-I2
C with KT
C KT = Triangle index and number of currently stored
C triangles
C LP = LIST pointer
C LP2 = Pointer to N2 as a neighbor of N1
C LPL = Pointer to the last neighbor of I1
C LPLN1 = Pointer to the last neighbor of N1
C N1,N2,N3 = Nodal indexes of triangle KT
C NM2 = N-2
C
C
C Test for invalid input parameters.
C
IF (N .LT. 3 .OR. (NROW .NE. 6 .AND. NROW .NE. 9))
. GO TO 11
C
C Initialize parameters for loop on triangles KT = (N1,N2,
C N3), where N1 < N2 and N1 < N3.
C
C ARCS = TRUE iff arc indexes are to be stored.
C KA,KT = Numbers of currently stored arcs and triangles.
C NM2 = Upper bound on candidates for N1.
C
ARCS = NROW .EQ. 9
KA = 0
KT = 0
NM2 = N-2
C
C Loop on nodes N1.
C
DO 9 N1 = 1,NM2
C
C Loop on pairs of adjacent neighbors (N2,N3). LPLN1 points
C to the last neighbor of N1, and LP2 points to N2.
C
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
C
C Add a new triangle KT = (N1,N2,N3).
C
KT = KT + 1
LTRI(1,KT) = N1
LTRI(2,KT) = N2
LTRI(3,KT) = N3
C
C Loop on triangle sides (I2,I1) with neighboring triangles
C KN = (I1,I2,I3).
C
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
C
C Set I3 to the neighbor of I1 that follows I2 unless
C I2->I1 is a boundary arc.
C
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
C
C I2 is the last neighbor of I1 unless the data structure
C is invalid. Bypass the search for a neighboring
C triangle if I2->I1 is a boundary arc.
C
IF (ABS(LIST(LP)) .NE. I2) GO TO 12
KN = 0
IF (LIST(LP) .LT. 0) GO TO 6
C
C I2->I1 is not a boundary arc, and LP points to I2 as
C a neighbor of I1.
C
3 LP = LPTR(LP)
I3 = ABS(LIST(LP))
C
C Find J such that LTRI(J,KN) = I3 (not used if KN > KT),
C and permute the vertex indexes of KN so that I1 is
C smallest.
C
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
C
C Test for KN > KT (triangle index not yet assigned).
C
IF (I1 .GT. N1) GO TO 7
C
C Find KN, if it exists, by searching the triangle list in
C reverse order.
C
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
C
C Store KT as a neighbor of KN.
C
5 LTRI(J+3,KN) = KT
C
C Store KN as a neighbor of KT, and add a new arc KA.
C
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
C
C Bottom of loop on triangles.
C
8 IF (LP2 .NE. LPLN1) GO TO 1
9 CONTINUE
C
C No errors encountered.
C
NT = KT
IER = 0
RETURN
C
C Invalid input parameter.
C
11 NT = 0
IER = 1
RETURN
C
C Invalid triangulation data structure: I1 is a neighbor of
C I2, but I2 is not a neighbor of I1.
C
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)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/02/98
C
C This subroutine prints the triangle list created by Sub-
C routine TRLIST and, optionally, the nodal coordinates
C (either latitude and longitude or Cartesian coordinates)
C on logical unit LOUT. The numbers of boundary nodes,
C triangles, and arcs are also printed.
C
C
C On input:
C
C N = Number of nodes in the triangulation.
C 3 .LE. N .LE. 9999.
C
C X,Y,Z = Arrays of length N containing the Cartesian
C coordinates of the nodes if IFLAG = 0, or
C (X and Y only) arrays of length N containing
C longitude and latitude, respectively, if
C IFLAG > 0, or unused dummy parameters if
C IFLAG < 0.
C
C IFLAG = Nodal coordinate option indicator:
C IFLAG = 0 if X, Y, and Z (assumed to contain
C Cartesian coordinates) are to be
C printed (to 6 decimal places).
C IFLAG > 0 if only X and Y (assumed to con-
C tain longitude and latitude) are
C to be printed (to 6 decimal
C places).
C IFLAG < 0 if only the adjacency lists are to
C be printed.
C
C NROW = Number of rows (entries per triangle) re-
C served for the triangle list LTRI. The value
C must be 6 if only the vertex indexes and
C neighboring triangle indexes are stored, or 9
C if arc indexes are also stored.
C
C NT = Number of triangles in the triangulation.
C 1 .LE. NT .LE. 9999.
C
C LTRI = NROW by NT array whose J-th column contains
C the vertex nodal indexes (first three rows),
C neighboring triangle indexes (second three
C rows), and, if NROW = 9, arc indexes (last
C three rows) associated with triangle J for
C J = 1,...,NT.
C
C LOUT = Logical unit number for output. If LOUT is
C not in the range 0 to 99, output is written
C to unit 6.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C The triangle list and nodal coordinates (as specified by
C IFLAG) are written to unit LOUT.
C
C Modules required by TRLPRT: None
C
C***********************************************************
C
INTEGER I, K, LUN, NA, NB, NL, NLMAX, NMAX
DATA NMAX/9999/, NLMAX/58/
C
C Local parameters:
C
C I = DO-loop, nodal index, and row index for LTRI
C K = DO-loop and triangle index
C LUN = Logical unit number for output
C NA = Number of triangulation arcs
C NB = Number of boundary nodes
C NL = Number of lines printed on the current page
C NLMAX = Maximum number of print lines per page (except
C for the last page which may have two addi-
C tional lines)
C NMAX = Maximum value of N and NT (4-digit format)
C
LUN = LOUT
IF (LUN .LT. 0 .OR. LUN .GT. 99) LUN = 6
C
C Print a heading and test for invalid input.
C
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
C
C Print an error message and exit.
C
WRITE (LUN,110) N, NROW, NT
RETURN
ENDIF
IF (IFLAG .EQ. 0) THEN
C
C Print X, Y, and Z.
C
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
C
C Print X (longitude) and Y (latitude).
C
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
C
C Print the triangulation LTRI.
C
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
C
C Print NB, NA, and NT (boundary nodes, arcs, and
C triangles).
C
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
C
C Print formats:
C
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)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/08/99
C
C This subroutine creates a Delaunay triangulation of a
C set of N arbitrarily distributed points, referred to as
C nodes, on the surface of the unit sphere. The Delaunay
C triangulation is defined as a set of (spherical) triangles
C with the following five properties:
C
C 1) The triangle vertices are nodes.
C 2) No triangle contains a node other than its vertices.
C 3) The interiors of the triangles are pairwise disjoint.
C 4) The union of triangles is the convex hull of the set
C of nodes (the smallest convex set that contains
C the nodes). If the nodes are not contained in a
C single hemisphere, their convex hull is the en-
C tire sphere and there are no boundary nodes.
C Otherwise, there are at least three boundary nodes.
C 5) The interior of the circumcircle of each triangle
C contains no node.
C
C The first four properties define a triangulation, and the
C last property results in a triangulation which is as close
C as possible to equiangular in a certain sense and which is
C uniquely defined unless four or more nodes lie in a common
C plane. This property makes the triangulation well-suited
C for solving closest-point problems and for triangle-based
C interpolation.
C
C Provided the nodes are randomly ordered, the algorithm
C has expected time complexity O(N*log(N)) for most nodal
C distributions. Note, however, that the complexity may be
C as high as O(N**2) if, for example, the nodes are ordered
C on increasing latitude.
C
C Spherical coordinates (latitude and longitude) may be
C converted to Cartesian coordinates by Subroutine TRANS.
C
C The following is a list of the software package modules
C which a user may wish to call directly:
C
C ADDNOD - Updates the triangulation by appending a new
C node.
C
C AREAS - Returns the area of a spherical triangle.
C
C BNODES - Returns an array containing the indexes of the
C boundary nodes (if any) in counterclockwise
C order. Counts of boundary nodes, triangles,
C and arcs are also returned.
C
C CIRCUM - Returns the circumcenter of a spherical trian-
C gle.
C
C CRLIST - Returns the set of triangle circumcenters
C (Voronoi vertices) and circumradii associated
C with a triangulation.
C
C DELARC - Deletes a boundary arc from a triangulation.
C
C DELNOD - Updates the triangulation with a nodal deletion.
C
C EDGE - Forces an arbitrary pair of nodes to be connec-
C ted by an arc in the triangulation.
C
C GETNP - Determines the ordered sequence of L closest
C nodes to a given node, along with the associ-
C ated distances.
C
C INSIDE - Locates a point relative to a polygon on the
C surface of the sphere.
C
C INTRSC - Returns the point of intersection between a
C pair of great circle arcs.
C
C JRAND - Generates a uniformly distributed pseudo-random
C integer.
C
C LEFT - Locates a point relative to a great circle.
C
C NEARND - Returns the index of the nearest node to an
C arbitrary point, along with its squared
C distance.
C
C SCOORD - Converts a point from Cartesian coordinates to
C spherical coordinates.
C
C STORE - Forces a value to be stored in main memory so
C that the precision of floating point numbers
C in memory locations rather than registers is
C computed.
C
C TRANS - Transforms spherical coordinates into Cartesian
C coordinates on the unit sphere for input to
C Subroutine TRMESH.
C
C TRLIST - Converts the triangulation data structure to a
C triangle list more suitable for use in a fin-
C ite element code.
C
C TRLPRT - Prints the triangle list created by Subroutine
C TRLIST.
C
C TRMESH - Creates a Delaunay triangulation of a set of
C nodes.
C
C TRPLOT - Creates a level-2 Encapsulated Postscript (EPS)
C file containing a triangulation plot.
C
C TRPRNT - Prints the triangulation data structure and,
C optionally, the nodal coordinates.
C
C VRPLOT - Creates a level-2 Encapsulated Postscript (EPS)
C file containing a Voronoi diagram plot.
C
C
C On input:
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y,Z = Arrays of length N containing the Cartesian
C coordinates of distinct nodes. (X(K),Y(K),
C Z(K)) is referred to as node K, and K is re-
C ferred to as a nodal index. It is required
C that X(K)**2 + Y(K)**2 + Z(K)**2 = 1 for all
C K. The first three nodes must not be col-
C linear (lie on a common great circle).
C
C The above parameters are not altered by this routine.
C
C LIST,LPTR = Arrays of length at least 6N-12.
C
C LEND = Array of length at least N.
C
C NEAR,NEXT,DIST = Work space arrays of length at
C least N. The space is used to
C efficiently determine the nearest
C triangulation node to each un-
C processed node for use by ADDNOD.
C
C On output:
C
C LIST = Set of nodal indexes which, along with LPTR,
C LEND, and LNEW, define the triangulation as a
C set of N adjacency lists -- counterclockwise-
C ordered sequences of neighboring nodes such
C that the first and last neighbors of a bound-
C ary node are boundary nodes (the first neigh-
C bor of an interior node is arbitrary). In
C order to distinguish between interior and
C boundary nodes, the last neighbor of each
C boundary node is represented by the negative
C of its index.
C
C LPTR = Set of pointers (LIST indexes) in one-to-one
C correspondence with the elements of LIST.
C LIST(LPTR(I)) indexes the node which follows
C LIST(I) in cyclical counterclockwise order
C (the first neighbor follows the last neigh-
C bor).
C
C LEND = Set of pointers to adjacency lists. LEND(K)
C points to the last neighbor of node K for
C K = 1,...,N. Thus, LIST(LEND(K)) < 0 if and
C only if K is a boundary node.
C
C LNEW = Pointer to the first empty location in LIST
C and LPTR (list length plus one). LIST, LPTR,
C LEND, and LNEW are not altered if IER < 0,
C and are incomplete if IER > 0.
C
C NEAR,NEXT,DIST = Garbage.
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = -1 if N < 3 on input.
C IER = -2 if the first three nodes are
C collinear.
C IER = L if nodes L and M coincide for some
C M > L. The data structure represents
C a triangulation of nodes 1 to M-1 in
C this case.
C
C Modules required by TRMESH: ADDNOD, BDYADD, COVSPH,
C INSERT, INTADD, JRAND,
C LEFT, LSTPTR, STORE, SWAP,
C SWPTST, TRFIND
C
C Intrinsic function called by TRMESH: ABS
C
C***********************************************************
C
INTEGER I, I0, J, K, LP, LPL, NEXTI, NN
LOGICAL LEFT
REAL D, D1, D2, D3
C
C Local parameters:
C
C D = (Negative cosine of) distance from node K to
C node I
C D1,D2,D3 = Distances from node K to nodes 1, 2, and 3,
C respectively
C I,J = Nodal indexes
C I0 = Index of the node preceding I in a sequence of
C unprocessed nodes: I = NEXT(I0)
C K = Index of node to be added and DO-loop index:
C K > 3
C LP = LIST index (pointer) of a neighbor of K
C LPL = Pointer to the last neighbor of K
C NEXTI = NEXT(I)
C NN = Local copy of N
C
NN = N
IF (NN .LT. 3) THEN
IER = -1
RETURN
ENDIF
C
C Store the first triangle in the linked list.
C
IF ( .NOT. LEFT (X(1),Y(1),Z(1),X(2),Y(2),Z(2),
. X(3),Y(3),Z(3)) ) THEN
C
C The first triangle is (3,2,1) = (2,1,3) = (1,3,2).
C
LIST(1) = 3
LPTR(1) = 2
LIST(2) = -2
LPTR(2) = 1
LEND(1) = 2
C
LIST(3) = 1
LPTR(3) = 4
LIST(4) = -3
LPTR(4) = 3
LEND(2) = 4
C
LIST(5) = 2
LPTR(5) = 6
LIST(6) = -1
LPTR(6) = 5
LEND(3) = 6
C
ELSEIF ( .NOT. LEFT(X(2),Y(2),Z(2),X(1),Y(1),Z(1),
. X(3),Y(3),Z(3)) )
. THEN
C
C The first triangle is (1,2,3): 3 Strictly Left 1->2,
C i.e., node 3 lies in the left hemisphere defined by
C arc 1->2.
C
LIST(1) = 2
LPTR(1) = 2
LIST(2) = -3
LPTR(2) = 1
LEND(1) = 2
C
LIST(3) = 3
LPTR(3) = 4
LIST(4) = -1
LPTR(4) = 3
LEND(2) = 4
C
LIST(5) = 1
LPTR(5) = 6
LIST(6) = -2
LPTR(6) = 5
LEND(3) = 6
C
ELSE
C
C The first three nodes are collinear.
C
IER = -2
RETURN
ENDIF
C
C Initialize LNEW and test for N = 3.
C
LNEW = 7
IF (NN .EQ. 3) THEN
IER = 0
RETURN
ENDIF
C
C A nearest-node data structure (NEAR, NEXT, and DIST) is
C used to obtain an expected-time (N*log(N)) incremental
C algorithm by enabling constant search time for locating
C each new node in the triangulation.
C
C For each unprocessed node K, NEAR(K) is the index of the
C triangulation node closest to K (used as the starting
C point for the search in Subroutine TRFIND) and DIST(K)
C is an increasing function of the arc length (angular
C distance) between nodes K and NEAR(K): -Cos(a) for arc
C length a.
C
C Since it is necessary to efficiently find the subset of
C unprocessed nodes associated with each triangulation
C node J (those that have J as their NEAR entries), the
C subsets are stored in NEAR and NEXT as follows: for
C each node J in the triangulation, I = NEAR(J) is the
C first unprocessed node in J's set (with I = 0 if the
C set is empty), L = NEXT(I) (if I > 0) is the second,
C NEXT(L) (if L > 0) is the third, etc. The nodes in each
C set are initially ordered by increasing indexes (which
C maximizes efficiency) but that ordering is not main-
C tained as the data structure is updated.
C
C Initialize the data structure for the single triangle.
C
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
C
C Add the remaining nodes
C
DO 6 K = 4,NN
CALL ADDNOD (NEAR(K),K,X,Y,Z, LIST,LPTR,LEND,
. LNEW, IER)
IF (IER .NE. 0) RETURN
C
C Remove K from the set of unprocessed nodes associated
C with NEAR(K).
C
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
C
C Loop on neighbors J of node K.
C
LPL = LEND(K)
LP = LPL
3 LP = LPTR(LP)
J = ABS(LIST(LP))
C
C Loop on elements I in the sequence of unprocessed nodes
C associated with J: K is a candidate for replacing J
C as the nearest triangulation node to I. The next value
C of I in the sequence, NEXT(I), must be saved before I
C is moved because it is altered by adding I to K's set.
C
I = NEAR(J)
4 IF (I .EQ. 0) GO TO 5
NEXTI = NEXT(I)
C
C Test for the distance from I to K less than the distance
C from I to J.
C
D = -(X(I)*X(K) + Y(I)*Y(K) + Z(I)*Z(K))
IF (D .LT. DIST(I)) THEN
C
C Replace J by K as the nearest triangulation node to I:
C update NEAR(I) and DIST(I), and remove I from J's set
C of unprocessed nodes and add it to K's set.
C
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
C
C Bottom of loop on I.
C
I = NEXTI
GO TO 4
C
C Bottom of loop on neighbors J.
C
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)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/16/98
C
C This subroutine creates a level-2 Encapsulated Post-
C script (EPS) file containing a graphical display of a
C triangulation of a set of nodes on the unit sphere. The
C visible nodes are projected onto the plane that contains
C the origin and has normal defined by a user-specified eye-
C position. Projections of adjacent (visible) nodes are
C connected by line segments.
C
C
C On input:
C
C LUN = Logical unit number in the range 0 to 99.
C The unit should be opened with an appropriate
C file name before the call to this routine.
C
C PLTSIZ = Plot size in inches. A circular window in
C the projection plane is mapped to a circu-
C lar viewport with diameter equal to .88*
C PLTSIZ (leaving room for labels outside the
C viewport). The viewport is centered on the
C 8.5 by 11 inch page, and its boundary is
C drawn. 1.0 .LE. PLTSIZ .LE. 8.5.
C
C ELAT,ELON = Latitude and longitude (in degrees) of
C the center of projection E (the center
C of the plot). The projection plane is
C the plane that contains the origin and
C has E as unit normal. In a rotated
C coordinate system for which E is the
C north pole, the projection plane con-
C tains the equator, and only northern
C hemisphere nodes are visible (from the
C point at infinity in the direction E).
C These are projected orthogonally onto
C the projection plane (by zeroing the z-
C component in the rotated coordinate
C system). ELAT and ELON must be in the
C range -90 to 90 and -180 to 180, respec-
C tively.
C
C A = Angular distance in degrees from E to the boun-
C dary of a circular window against which the
C triangulation is clipped. The projected window
C is a disk of radius r = Sin(A) centered at the
C origin, and only visible nodes whose projections
C are within distance r of the origin are included
C in the plot. Thus, if A = 90, the plot includes
C the entire hemisphere centered at E. 0 .LT. A
C .LE. 90.
C
C N = Number of nodes in the triangulation. N .GE. 3.
C
C X,Y,Z = Arrays of length N containing the Cartesian
C coordinates of the nodes (unit vectors).
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C TITLE = Type CHARACTER variable or constant contain-
C ing a string to be centered above the plot.
C The string must be enclosed in parentheses;
C i.e., the first and last characters must be
C '(' and ')', respectively, but these are not
C displayed. TITLE may have at most 80 char-
C acters including the parentheses.
C
C NUMBR = Option indicator: If NUMBR = TRUE, the
C nodal indexes are plotted next to the nodes.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if LUN, PLTSIZ, or N is outside its
C valid range.
C IER = 2 if ELAT, ELON, or A is outside its
C valid range.
C IER = 3 if an error was encountered in writing
C to unit LUN.
C
C The values in the data statement below may be altered
C in order to modify various plotting options.
C
C Modules required by TRPLOT: None
C
C Intrinsic functions called by TRPLOT: ABS, ATAN, COS,
C NINT, REAL, SIN,
C SQRT
C
C***********************************************************
C
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
C
DATA ANNOT/.TRUE./, FSIZN/10.0/, FSIZT/16.0/
C
C Local parameters:
C
C ANNOT = Logical variable with value TRUE iff the plot
C is to be annotated with the values of ELAT,
C ELON, and A
C CF = Conversion factor for degrees to radians
C CT = Cos(ELAT)
C EX,EY,EZ = Cartesian coordinates of the eye-position E
C FSIZN = Font size in points for labeling nodes with
C their indexes if NUMBR = TRUE
C FSIZT = Font size in points for the title (and
C annotation if ANNOT = TRUE)
C IPX1,IPY1 = X and y coordinates (in points) of the lower
C left corner of the bounding box or viewport
C box
C IPX2,IPY2 = X and y coordinates (in points) of the upper
C right corner of the bounding box or viewport
C box
C IR = Half the width (height) of the bounding box or
C viewport box in points -- viewport radius
C LP = LIST index (pointer)
C LPL = Pointer to the last neighbor of N0
C N0 = Index of a node whose incident arcs are to be
C drawn
C N1 = Neighbor of N0
C R11...R23 = Components of the first two rows of a rotation
C that maps E to the north pole (0,0,1)
C SF = Scale factor for mapping world coordinates
C (window coordinates in [-WR,WR] X [-WR,WR])
C to viewport coordinates in [IPX1,IPX2] X
C [IPY1,IPY2]
C T = Temporary variable
C TX,TY = Translation vector for mapping world coordi-
C nates to viewport coordinates
C WR = Window radius r = Sin(A)
C WRS = WR**2
C X0,Y0,Z0 = Coordinates of N0 in the rotated coordinate
C system or label location (X0,Y0)
C X1,Y1,Z1 = Coordinates of N1 in the rotated coordinate
C system or intersection of edge N0-N1 with
C the equator (in the rotated coordinate
C system)
C
C
C Test for invalid parameters.
C
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
C
C Compute a conversion factor CF for degrees to radians
C and compute the window radius WR.
C
CF = ATAN(1.0)/45.0
WR = SIN(CF*A)
WRS = WR*WR
C
C Compute the lower left (IPX1,IPY1) and upper right
C (IPX2,IPY2) corner coordinates of the bounding box.
C The coordinates, specified in default user space units
C (points, at 72 points/inch with origin at the lower
C left corner of the page), are chosen to preserve the
C square aspect ratio, and to center the plot on the 8.5
C by 11 inch page. The center of the page is (306,396),
C and IR = PLTSIZ/2 in points.
C
IR = NINT(36.0*PLTSIZ)
IPX1 = 306 - IR
IPX2 = 306 + IR
IPY1 = 396 - IR
IPY2 = 396 + IR
C
C Output header comments.
C
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')
C
C Set (IPX1,IPY1) and (IPX2,IPY2) to the corner coordinates
C of a viewport box obtained by shrinking the bounding box
C by 12% in each dimension.
C
IR = NINT(0.88*REAL(IR))
IPX1 = 306 - IR
IPX2 = 306 + IR
IPY1 = 396 - IR
IPY2 = 396 + IR
C
C Set the line thickness to 2 points, and draw the
C viewport boundary.
C
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')
C
C Set up an affine mapping from the window box [-WR,WR] X
C [-WR,WR] to the viewport box.
C
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')
C
C The line thickness must be changed to reflect the new
C scaling which is applied to all subsequent output.
C Set it to 1.0 point.
C
T = 1.0/SF
WRITE (LUN,110,ERR=13) T
C
C Save the current graphics state, and set the clip path to
C the boundary of the window.
C
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')
C
C Compute the Cartesian coordinates of E and the components
C of a rotation R which maps E to the north pole (0,0,1).
C R is taken to be a rotation about the z-axis (into the
C yz-plane) followed by a rotation about the x-axis chosen
C so that the view-up direction is (0,0,1), or (-1,0,0) if
C E is the north or south pole.
C
C ( R11 R12 0 )
C R = ( R21 R22 R23 )
C ( EX EY EZ )
C
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
C
C Loop on visible nodes N0 that project to points (X0,Y0) in
C the window.
C
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
C
C Loop on neighbors N1 of N0. LPL points to the last
C neighbor of N0. Copy the components of N1 into P.
C
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
C
C N1 is a 'southern hemisphere' point. Move it to the
C intersection of edge N0-N1 with the equator so that
C the edge is clipped properly. Z1 is implicitly set
C to 0.
C
X1 = Z0*X1 - Z1*X0
Y1 = Z0*Y1 - Z1*Y0
T = SQRT(X1*X1+Y1*Y1)
X1 = X1/T
Y1 = Y1/T
ENDIF
C
C If node N1 is in the window and N1 < N0, bypass edge
C N0->N1 (since edge N1->N0 has already been drawn).
C
IF ( Z1 .GE. 0.0 .AND. X1*X1 + Y1*Y1 .LE. WRS
. .AND. N1 .LT. N0 ) GO TO 2
C
C Add the edge to the path.
C
WRITE (LUN,180,ERR=13) X0, Y0, X1, Y1
180 FORMAT (2F12.6,' moveto',2F12.6,' lineto')
C
C Bottom of loops.
C
2 IF (LP .NE. LPL) GO TO 1
3 CONTINUE
C
C Paint the path and restore the saved graphics state (with
C no clip path).
C
WRITE (LUN,130,ERR=13)
WRITE (LUN,190,ERR=13)
190 FORMAT ('grestore')
IF (NUMBR) THEN
C
C Nodes in the window are to be labeled with their indexes.
C Convert FSIZN from points to world coordinates, and
C output the commands to select a font and scale it.
C
T = FSIZN/SF
WRITE (LUN,200,ERR=13) T
200 FORMAT ('/Helvetica findfont'/
. F12.6,' scalefont setfont')
C
C Loop on visible nodes N0 that project to points (X0,Y0) in
C the window.
C
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
C
C Move to (X0,Y0) and draw the label N0. The first char-
C acter will will have its lower left corner about one
C character width to the right of the nodal position.
C
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
C
C Convert FSIZT from points to world coordinates, and output
C the commands to select a font and scale it.
C
T = FSIZT/SF
WRITE (LUN,200,ERR=13) T
C
C Display TITLE centered above the plot:
C
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
C
C Display the window center and radius below the plot.
C
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
C
C Paint the path and output the showpage command and
C end-of-file indicator.
C
WRITE (LUN,270,ERR=13)
270 FORMAT ('stroke'/
. 'showpage'/
. '%%EOF')
C
C HP's interpreters require a one-byte End-of-PostScript-Job
C indicator (to eliminate a timeout error message):
C ASCII 4.
C
WRITE (LUN,280,ERR=13) CHAR(4)
280 FORMAT (A1)
C
C No error encountered.
C
IER = 0
RETURN
C
C Invalid input parameter LUN, PLTSIZ, or N.
C
11 IER = 1
RETURN
C
C Invalid input parameter ELAT, ELON, or A.
C
12 IER = 2
RETURN
C
C Error writing to unit LUN.
C
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)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/25/98
C
C This subroutine prints the triangulation adjacency lists
C created by Subroutine TRMESH and, optionally, the nodal
C coordinates (either latitude and longitude or Cartesian
C coordinates) on logical unit LOUT. The list of neighbors
C of a boundary node is followed by index 0. The numbers of
C boundary nodes, triangles, and arcs are also printed.
C
C
C On input:
C
C N = Number of nodes in the triangulation. N .GE. 3
C and N .LE. 9999.
C
C X,Y,Z = Arrays of length N containing the Cartesian
C coordinates of the nodes if IFLAG = 0, or
C (X and Y only) arrays of length N containing
C longitude and latitude, respectively, if
C IFLAG > 0, or unused dummy parameters if
C IFLAG < 0.
C
C IFLAG = Nodal coordinate option indicator:
C IFLAG = 0 if X, Y, and Z (assumed to contain
C Cartesian coordinates) are to be
C printed (to 6 decimal places).
C IFLAG > 0 if only X and Y (assumed to con-
C tain longitude and latitude) are
C to be printed (to 6 decimal
C places).
C IFLAG < 0 if only the adjacency lists are to
C be printed.
C
C LIST,LPTR,LEND = Data structure defining the trian-
C gulation. Refer to Subroutine
C TRMESH.
C
C LOUT = Logical unit for output. If LOUT is not in
C the range 0 to 99, output is written to
C logical unit 6.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C The adjacency lists and nodal coordinates (as specified
C by IFLAG) are written to unit LOUT.
C
C Modules required by TRPRNT: None
C
C***********************************************************
C
INTEGER I, INC, K, LP, LPL, LUN, NA, NABOR(400), NB,
. ND, NL, NLMAX, NMAX, NODE, NN, NT
DATA NMAX/9999/, NLMAX/58/
C
C Local parameters:
C
C I = NABOR index (1 to K)
C INC = Increment for NL associated with an adjacency list
C K = Counter and number of neighbors of NODE
C LP = LIST pointer of a neighbor of NODE
C LPL = Pointer to the last neighbor of NODE
C LUN = Logical unit for output (copy of LOUT)
C NA = Number of arcs in the triangulation
C NABOR = Array containing the adjacency list associated
C with NODE, with zero appended if NODE is a
C boundary node
C NB = Number of boundary nodes encountered
C ND = Index of a neighbor of NODE (or negative index)
C NL = Number of lines that have been printed on the
C current page
C NLMAX = Maximum number of print lines per page (except
C for the last page which may have two addi-
C tional lines)
C NMAX = Upper bound on N (allows 4-digit indexes)
C NODE = Index of a node and DO-loop index (1 to N)
C NN = Local copy of N
C NT = Number of triangles in the triangulation
C
NN = N
LUN = LOUT
IF (LUN .LT. 0 .OR. LUN .GT. 99) LUN = 6
C
C Print a heading and test the range of N.
C
WRITE (LUN,100) NN
IF (NN .LT. 3 .OR. NN .GT. NMAX) THEN
C
C N is outside its valid range.
C
WRITE (LUN,110)
RETURN
ENDIF
C
C Initialize NL (the number of lines printed on the current
C page) and NB (the number of boundary nodes encountered).
C
NL = 6
NB = 0
IF (IFLAG .LT. 0) THEN
C
C Print LIST only. K is the number of neighbors of NODE
C that have been stored in NABOR.
C
WRITE (LUN,101)
DO 2 NODE = 1,NN
LPL = LEND(NODE)
LP = LPL
K = 0
C
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
C
C NODE is a boundary node. Correct the sign of the last
C neighbor, add 0 to the end of the list, and increment
C NB.
C
NABOR(K) = -ND
K = K + 1
NABOR(K) = 0
NB = NB + 1
ENDIF
C
C Increment NL and print the list of neighbors.
C
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
C
C Print X (longitude), Y (latitude), and LIST.
C
WRITE (LUN,102)
DO 4 NODE = 1,NN
LPL = LEND(NODE)
LP = LPL
K = 0
C
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
C
C NODE is a boundary node.
C
NABOR(K) = -ND
K = K + 1
NABOR(K) = 0
NB = NB + 1
ENDIF
C
C Increment NL and print X, Y, and NABOR.
C
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
C
C Print X, Y, Z, and LIST.
C
WRITE (LUN,103)
DO 6 NODE = 1,NN
LPL = LEND(NODE)
LP = LPL
K = 0
C
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
C
C NODE is a boundary node.
C
NABOR(K) = -ND
K = K + 1
NABOR(K) = 0
NB = NB + 1
ENDIF
C
C Increment NL and print X, Y, Z, and NABOR.
C
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
C
C Print NB, NA, and NT (boundary nodes, arcs, and
C triangles).
C
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
C
C Print formats:
C
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)
C
C***********************************************************
C
C From STRIPACK
C Robert J. Renka
C Dept. of Computer Science
C Univ. of North Texas
C renka@cs.unt.edu
C 07/16/98
C
C This subroutine creates a level-2 Encapsulated Post-
C script (EPS) file containing a graphical depiction of a
C Voronoi diagram of a set of nodes on the unit sphere.
C The visible vertices are projected onto the plane that
C contains the origin and has normal defined by a user-
C specified eye-position. Projections of adjacent (visible)
C Voronoi vertices are connected by line segments.
C
C The parameters defining the Voronoi diagram may be com-
C puted by Subroutine CRLIST.
C
C
C On input:
C
C LUN = Logical unit number in the range 0 to 99.
C The unit should be opened with an appropriate
C file name before the call to this routine.
C
C PLTSIZ = Plot size in inches. A circular window in
C the projection plane is mapped to a circu-
C lar viewport with diameter equal to .88*
C PLTSIZ (leaving room for labels outside the
C viewport). The viewport is centered on the
C 8.5 by 11 inch page, and its boundary is
C drawn. 1.0 .LE. PLTSIZ .LE. 8.5.
C
C ELAT,ELON = Latitude and longitude (in degrees) of
C the center of projection E (the center
C of the plot). The projection plane is
C the plane that contains the origin and
C has E as unit normal. In a rotated
C coordinate system for which E is the
C north pole, the projection plane con-
C tains the equator, and only northern
C hemisphere points are visible (from the
C point at infinity in the direction E).
C These are projected orthogonally onto
C the projection plane (by zeroing the z-
C component in the rotated coordinate
C system). ELAT and ELON must be in the
C range -90 to 90 and -180 to 180, respec-
C tively.
C
C A = Angular distance in degrees from E to the boun-
C dary of a circular window against which the
C Voronoi diagram is clipped. The projected win-
C dow is a disk of radius r = Sin(A) centered at
C the origin, and only visible vertices whose
C projections are within distance r of the origin
C are included in the plot. Thus, if A = 90, the
C plot includes the entire hemisphere centered at
C E. 0 .LT. A .LE. 90.
C
C N = Number of nodes (Voronoi centers) and Voronoi
C regions. N .GE. 3.
C
C X,Y,Z = Arrays of length N containing the Cartesian
C coordinates of the nodes (unit vectors).
C
C NT = Number of Voronoi region vertices (triangles,
C including those in the extended triangulation
C if the number of boundary nodes NB is nonzero):
C NT = 2*N-4.
C
C LISTC = Array of length 3*NT containing triangle
C indexes (indexes to XC, YC, and ZC) stored
C in 1-1 correspondence with LIST/LPTR entries
C (or entries that would be stored in LIST for
C the extended triangulation): the index of
C triangle (N1,N2,N3) is stored in LISTC(K),
C LISTC(L), and LISTC(M), where LIST(K),
C LIST(L), and LIST(M) are the indexes of N2
C as a neighbor of N1, N3 as a neighbor of N2,
C and N1 as a neighbor of N3. The Voronoi
C region associated with a node is defined by
C the CCW-ordered sequence of circumcenters in
C one-to-one correspondence with its adjacency
C list (in the extended triangulation).
C
C LPTR = Array of length 3*NT = 6*N-12 containing a
C set of pointers (LISTC indexes) in one-to-one
C correspondence with the elements of LISTC.
C LISTC(LPTR(I)) indexes the triangle which
C follows LISTC(I) in cyclical counterclockwise
C order (the first neighbor follows the last
C neighbor).
C
C LEND = Array of length N containing a set of
C pointers to triangle lists. LP = LEND(K)
C points to a triangle (indexed by LISTC(LP))
C containing node K for K = 1 to N.
C
C XC,YC,ZC = Arrays of length NT containing the
C Cartesian coordinates of the triangle
C circumcenters (Voronoi vertices).
C XC(I)**2 + YC(I)**2 + ZC(I)**2 = 1.
C
C TITLE = Type CHARACTER variable or constant contain-
C ing a string to be centered above the plot.
C The string must be enclosed in parentheses;
C i.e., the first and last characters must be
C '(' and ')', respectively, but these are not
C displayed. TITLE may have at most 80 char-
C acters including the parentheses.
C
C NUMBR = Option indicator: If NUMBR = TRUE, the
C nodal indexes are plotted at the Voronoi
C region centers.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C IER = Error indicator:
C IER = 0 if no errors were encountered.
C IER = 1 if LUN, PLTSIZ, N, or NT is outside
C its valid range.
C IER = 2 if ELAT, ELON, or A is outside its
C valid range.
C IER = 3 if an error was encountered in writing
C to unit LUN.
C
C Modules required by VRPLOT: None
C
C Intrinsic functions called by VRPLOT: ABS, ATAN, COS,
C NINT, REAL, SIN,
C SQRT
C
C***********************************************************
C
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
C
DATA ANNOT/.TRUE./, FSIZN/10.0/, FSIZT/16.0/
C
C Local parameters:
C
C ANNOT = Logical variable with value TRUE iff the plot
C is to be annotated with the values of ELAT,
C ELON, and A
C CF = Conversion factor for degrees to radians
C CT = Cos(ELAT)
C EX,EY,EZ = Cartesian coordinates of the eye-position E
C FSIZN = Font size in points for labeling nodes with
C their indexes if NUMBR = TRUE
C FSIZT = Font size in points for the title (and
C annotation if ANNOT = TRUE)
C IN1,IN2 = Logical variables with value TRUE iff the
C projections of vertices KV1 and KV2, respec-
C tively, are inside the window
C IPX1,IPY1 = X and y coordinates (in points) of the lower
C left corner of the bounding box or viewport
C box
C IPX2,IPY2 = X and y coordinates (in points) of the upper
C right corner of the bounding box or viewport
C box
C IR = Half the width (height) of the bounding box or
C viewport box in points -- viewport radius
C KV1,KV2 = Endpoint indexes of a Voronoi edge
C LP = LIST index (pointer)
C LPL = Pointer to the last neighbor of N0
C N0 = Index of a node
C R11...R23 = Components of the first two rows of a rotation
C that maps E to the north pole (0,0,1)
C SF = Scale factor for mapping world coordinates
C (window coordinates in [-WR,WR] X [-WR,WR])
C to viewport coordinates in [IPX1,IPX2] X
C [IPY1,IPY2]
C T = Temporary variable
C TX,TY = Translation vector for mapping world coordi-
C nates to viewport coordinates
C WR = Window radius r = Sin(A)
C WRS = WR**2
C X0,Y0 = Projection plane coordinates of node N0 or
C label location
C X1,Y1,Z1 = Coordinates of vertex KV1 in the rotated
C coordinate system
C X2,Y2,Z2 = Coordinates of vertex KV2 in the rotated
C coordinate system or intersection of edge
C KV1-KV2 with the equator (in the rotated
C coordinate system)
C
C
C Test for invalid parameters.
C
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
C
C Compute a conversion factor CF for degrees to radians
C and compute the window radius WR.
C
CF = ATAN(1.0)/45.0
WR = SIN(CF*A)
WRS = WR*WR
C
C Compute the lower left (IPX1,IPY1) and upper right
C (IPX2,IPY2) corner coordinates of the bounding box.
C The coordinates, specified in default user space units
C (points, at 72 points/inch with origin at the lower
C left corner of the page), are chosen to preserve the
C square aspect ratio, and to center the plot on the 8.5
C by 11 inch page. The center of the page is (306,396),
C and IR = PLTSIZ/2 in points.
C
IR = NINT(36.0*PLTSIZ)
IPX1 = 306 - IR
IPX2 = 306 + IR
IPY1 = 396 - IR
IPY2 = 396 + IR
C
C Output header comments.
C
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')
C
C Set (IPX1,IPY1) and (IPX2,IPY2) to the corner coordinates
C of a viewport box obtained by shrinking the bounding box
C by 12% in each dimension.
C
IR = NINT(0.88*REAL(IR))
IPX1 = 306 - IR
IPX2 = 306 + IR
IPY1 = 396 - IR
IPY2 = 396 + IR
C
C Set the line thickness to 2 points, and draw the
C viewport boundary.
C
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')
C
C Set up an affine mapping from the window box [-WR,WR] X
C [-WR,WR] to the viewport box.
C
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')
C
C The line thickness must be changed to reflect the new
C scaling which is applied to all subsequent output.
C Set it to 1.0 point.
C
T = 1.0/SF
WRITE (LUN,110,ERR=13) T
C
C Save the current graphics state, and set the clip path to
C the boundary of the window.
C
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')
C
C Compute the Cartesian coordinates of E and the components
C of a rotation R which maps E to the north pole (0,0,1).
C R is taken to be a rotation about the z-axis (into the
C yz-plane) followed by a rotation about the x-axis chosen
C so that the view-up direction is (0,0,1), or (-1,0,0) if
C E is the north or south pole.
C
C ( R11 R12 0 )
C R = ( R21 R22 R23 )
C ( EX EY EZ )
C
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
C
C Loop on nodes (Voronoi centers) N0.
C LPL indexes the last neighbor of N0.
C
DO 3 N0 = 1,N
LPL = LEND(N0)
C
C Set KV2 to the first (and last) vertex index and compute
C its coordinates (X2,Y2,Z2) in the rotated coordinate
C system.
C
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)
C
C IN2 = TRUE iff KV2 is in the window.
C
IN2 = Z2 .GE. 0. .AND. X2*X2 + Y2*Y2 .LE. WRS
C
C Loop on neighbors N1 of N0. For each triangulation edge
C N0-N1, KV1-KV2 is the corresponding Voronoi edge.
C
LP = LPL
1 LP = LPTR(LP)
KV1 = KV2
X1 = X2
Y1 = Y2
Z1 = Z2
IN1 = IN2
KV2 = LISTC(LP)
C
C Compute the new values of (X2,Y2,Z2) and IN2.
C
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
C
C Add edge KV1-KV2 to the path iff both endpoints are inside
C the window and KV2 > KV1, or KV1 is inside and KV2 is
C outside (so that the edge is drawn only once).
C
IF (.NOT. IN1 .OR. (IN2 .AND. KV2 .LE. KV1))
. GO TO 2
IF (Z2 .LT. 0.) THEN
C
C KV2 is a 'southern hemisphere' point. Move it to the
C intersection of edge KV1-KV2 with the equator so that
C the edge is clipped properly. Z2 is implicitly set
C to 0.
C
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')
C
C Bottom of loops.
C
2 IF (LP .NE. LPL) GO TO 1
3 CONTINUE
C
C Paint the path and restore the saved graphics state (with
C no clip path).
C
WRITE (LUN,130,ERR=13)
WRITE (LUN,190,ERR=13)
190 FORMAT ('grestore')
IF (NUMBR) THEN
C
C Nodes in the window are to be labeled with their indexes.
C Convert FSIZN from points to world coordinates, and
C output the commands to select a font and scale it.
C
T = FSIZN/SF
WRITE (LUN,200,ERR=13) T
200 FORMAT ('/Helvetica findfont'/
. F12.6,' scalefont setfont')
C
C Loop on visible nodes N0 that project to points (X0,Y0) in
C the window.
C
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
C
C Move to (X0,Y0), and draw the label N0 with the origin
C of the first character at (X0,Y0).
C
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
C
C Convert FSIZT from points to world coordinates, and output
C the commands to select a font and scale it.
C
T = FSIZT/SF
WRITE (LUN,200,ERR=13) T
C
C Display TITLE centered above the plot:
C
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
C
C Display the window center and radius below the plot.
C
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
C
C Paint the path and output the showpage command and
C end-of-file indicator.
C
WRITE (LUN,270,ERR=13)
270 FORMAT ('stroke'/
. 'showpage'/
. '%%EOF')
C
C HP's interpreters require a one-byte End-of-PostScript-Job
C indicator (to eliminate a timeout error message):
C ASCII 4.
C
WRITE (LUN,280,ERR=13) CHAR(4)
280 FORMAT (A1)
C
C No error encountered.
C
IER = 0
RETURN
C
C Invalid input parameter LUN, PLTSIZ, N, or NT.
C
11 IER = 1
RETURN
C
C Invalid input parameter ELAT, ELON, or A.
C
12 IER = 2
RETURN
C
C Error writing to unit LUN.
C
13 IER = 3
RETURN
END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
# End of shell archive
exit 0