C ALGORITHM 623 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 4, C DEC., 1984, P. 437. SUBROUTINE ADNODE (KK,X,Y,Z, IADJ,IEND, IER) INTEGER KK, IADJ(1), IEND(KK), IER REAL X(KK), Y(KK), Z(KK) LOGICAL SWPTST EXTERNAL INDEX C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE ADDS NODE KK TO A TRIANGULATION OF THE C CONVEX HULL OF NODES 1,...,KK-1, PRODUCING A TRIANGULATION C OF THE CONVEX HULL OF NODES 1,...,KK. A SEQUENCE OF EDGE C SWAPS IS THEN APPLIED TO THE MESH, RESULTING IN AN OPTIMAL C TRIANGULATION. ADNODE IS PART OF AN INTERPOLATION PACKAGE C WHICH ALSO PROVIDES ROUTINES TO INITIALIZE THE DATA STRUC- C TURE, PLOT THE MESH, AND DELETE ARCS. C C INPUT PARAMETERS - KK - INDEX OF THE NODE TO BE ADDED C TO THE MESH. KK .GE. 4. C C X,Y,Z - VECTORS OF LENGTH .GE. KK CON- C TAINING CARTESIAN COORDINATES C OF THE NODES. (X(I),Y(I),Z(I)) C DEFINES NODE I FOR I = 1,...,KK. C C IADJ - SET OF ADJACENCY LISTS OF NODES C 1,...,KK-1. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IADJ AND IEND MAY BE CREATED BY TRMESH. C C KK, X, Y, AND Z ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION C OF NODE KK AS THE LAST C ENTRY. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS C WERE ENCOUNTERED. C IER = 1 IF ALL NODES C (INCLUDING KK) ARE C COLLINEAR. C C MODULES REFERENCED BY ADNODE - TRFIND, INTADD, BDYADD, C COVSPH, SHIFTD, INDEX, C SWPTST, SWAP C C*********************************************************** C INTEGER K, KM1, I1, I2, I3, INDKF, INDKL, NABOR1, . IO1, IO2, IN1, INDK1, IND2F, IND21 REAL P(3), DUM C C LOCAL PARAMETERS - C C K = LOCAL COPY OF KK C KM1 = K - 1 C I1,I2,I3 = VERTICES OF A TRIANGLE CONTAINING K C INDKF = IADJ INDEX OF THE FIRST NEIGHBOR OF K C INDKL = IADJ INDEX OF THE LAST NEIGHBOR OF K C NABOR1 = FIRST NEIGHBOR OF K BEFORE ANY SWAPS OCCUR C IO1,IO2 = ADJACENT NEIGHBORS OF K DEFINING AN ARC TO C BE TESTED FOR A SWAP C IN1 = VERTEX OPPOSITE K -- FIRST NEIGHBOR OF IO2 C WHICH PRECEDES IO1. IN1,IO1,IO2 ARE IN C COUNTERCLOCKWISE ORDER. C INDK1 = INDEX OF IO1 IN THE ADJACENCY LIST FOR K C IND2F = INDEX OF THE FIRST NEIGHBOR OF IO2 C IND21 = INDEX OF IO1 IN THE ADJACENCY LIST FOR IO2 C P = CARTESIAN COORDINATES OF NODE KK C DUM = DUMMY PARAMETER FOR CALL TO TRFIND C IER = 0 K = KK C C INITIALIZATION C KM1 = K - 1 P(1) = X(K) P(2) = Y(K) P(3) = Z(K) C C ADD NODE K TO THE MESH C CALL TRFIND(KM1,P,X,Y,Z,IADJ,IEND, DUM,DUM,DUM, . I1,I2,I3) IF (I1 .EQ. 0) GO TO 5 IF (I3 .NE. 0) CALL INTADD(K,I1,I2,I3, IADJ,IEND ) IF (I3 .EQ. 0 .AND. I1 .NE. I2) . CALL BDYADD (K,I1,I2, IADJ,IEND ) IF (I1 .EQ. I2) CALL COVSPH(K,I1, IADJ,IEND ) C C INITIALIZE VARIABLES FOR OPTIMIZATION OF THE MESH C INDKF = IEND(KM1) + 1 INDKL = IEND(K) NABOR1 = IADJ(INDKF) IO2 = NABOR1 INDK1 = INDKF + 1 IO1 = IADJ(INDK1) C C BEGIN LOOP -- FIND THE VERTEX OPPOSITE K C 1 IND2F = 1 IF (IO2 .NE. 1) IND2F = IEND(IO2-1) + 1 IND21 = INDEX(IO2,IO1,IADJ,IEND) IF (IND2F .EQ. IND21) GO TO 2 IN1 = IADJ(IND21-1) GO TO 3 C C IN1 IS THE LAST NEIGHBOR OF IO2 C 2 IND21 = IEND(IO2) IN1 = IADJ(IND21) IF (IN1 .EQ. 0) GO TO 4 C C SWAP TEST -- IF A SWAP OCCURS, TWO NEW ARCS ARE OPPOSITE K C AND MUST BE TESTED. INDK1 AND INDKF MUST BE C DECREMENTED. C 3 IF ( .NOT. SWPTST(IO1,IO2,IN1,K,X,Y,Z) ) GO TO 4 CALL SWAP(IN1,K,IO1,IO2, IADJ,IEND ) IO1 = IN1 INDK1 = INDK1 - 1 INDKF = INDKF - 1 GO TO 1 C C NO SWAP OCCURRED. RESET IO2 AND IO1, AND TEST FOR C TERMINATION. C 4 IF (IO1 .EQ. NABOR1) RETURN IO2 = IO1 INDK1 = INDK1 + 1 IF (INDK1 .GT. INDKL) INDK1 = INDKF IO1 = IADJ(INDK1) IF (IO1 .NE. 0) GO TO 1 RETURN C C ALL NODES ARE COLLINEAR C 5 IER = 1 RETURN END SUBROUTINE APLYR (X,Y,Z,CX,SX,CY,SY, XP,YP,ZP) REAL X, Y, Z, CX, SX, CY, SY, XP, YP, ZP C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE APPLIES THE ROTATION R DEFINED BY SUB- C ROUTINE CONSTR TO THE UNIT VECTOR (X Y Z)**T, I.E. (X,Y,Z) C IS ROTATED TO (XP,YP,ZP). IF (XP,YP,ZP) LIES IN THE C SOUTHERN HEMISPHERE (ZP .LT. 0), (XP,YP) ARE SET TO THE C COORDINATES OF THE NEAREST POINT OF THE EQUATOR, ZP RE- C MAINING UNCHANGED. C C INPUT PARAMETERS - X,Y,Z - COORDINATES OF A POINT ON THE C UNIT SPHERE. C C CX,SX,CY,SY - ELEMENTS OF THE ROTATION DE- C FINED BY CONSTR. C C INPUT PARAMETERS ARE NOT ALTERED EXCEPT AS NOTED BELOW. C C OUTPUT PARAMETERS - XP,YP,ZP - COORDINATES OF THE ROTATED C POINT ON THE SPHERE UNLESS C ZP .LT. 0, IN WHICH CASE C (XP,YP,0) IS THE CLOSEST C POINT OF THE EQUATOR TO THE C ROTATED POINT. STORAGE FOR C XP, YP, AND ZP MAY COINCIDE C WITH STORAGE FOR X, Y, AND C Z, RESPECTIVELY, IF THE C LATTER NEED NOT BE SAVED. C C MODULES REFERENCED BY APLYR - NONE C C INTRINSIC FUNCTION CALLED BY APLYR - SQRT C C*********************************************************** C REAL T C C LOCAL PARAMETER - C C T = TEMPORARY VARIABLE C T = SX*Y + CX*Z YP = CX*Y - SX*Z ZP = SY*X + CY*T XP = CY*X - SY*T IF (ZP .GE. 0.) RETURN C C MOVE (XP,YP,ZP) TO THE EQUATOR C T = SQRT(XP*XP + YP*YP) IF (T .EQ. 0.) GO TO 1 XP = XP/T YP = YP/T RETURN C C MOVE THE SOUTH POLE TO AN ARBITRARY POINT OF THE EQUATOR C 1 XP = 1. YP = 0. RETURN END SUBROUTINE APLYRT (G1P,G2P,CX,SX,CY,SY, G) REAL G1P, G2P, CX, SX, CY, SY, G(3) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE APPLIES THE INVERSE (TRANSPOSE) OF THE C ROTATION DEFINED BY SUBROUTINE CONSTR TO THE VECTOR C (G1P G2P 0)**T, I.E. THE GRADIENT (G1P,G2P,0) IN THE ROT- C ATED COORDINATE SYSTEM IS MAPPED TO (G1,G2,G3) IN THE C ORIGINAL COORDINATE SYSTEM. C C INPUT PARAMETERS - G1P,G2P - X- AND Y-COMPONENTS, RESPECT- C IVELY, OF THE GRADIENT IN THE C ROTATED COORDINATE SYSTEM. C C CX,SX,CY,SY - ELEMENTS OF THE ROTATION R C CONSTRUCTED BY SUBROUTINE C CONSTR. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - G - X-, Y-, AND Z-COMPONENTS (IN THAT C ORDER) OF THE INVERSE ROTATION C APPLIED TO (G1P,G2P,0) -- GRADIENT C IN THE ORIGINAL COORDINATE SYSTEM. C C MODULES REFERENCED BY APLYRT - NONE C C*********************************************************** C REAL T C C LOCAL PARAMETERS - C C T = TEMPORARY VARIABLE C T = SY*G1P G(1) = CY*G1P G(2) = CX*G2P - SX*T G(3) = -SX*G2P - CX*T RETURN END SUBROUTINE ARCINT (P,P1,P2,W1,W2,G1,G2, W,G,GN) REAL P(3), P1(3), P2(3), W1, W2, G1(3), G2(3), . W, G(3), GN C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN 3 POINTS P, P1, AND P2 LYING ON A COMMON GEODESIC C OF THE UNIT SPHERE WITH P BETWEEN P1 AND P2, ALONG WITH C DATA VALUES AND GRADIENTS AT P1 AND P2, THIS SUBROUTINE C COMPUTES AN INTERPOLATED VALUE W AND A GRADIENT VECTOR G C AT P. W IS COMPUTED BY HERMITE CUBIC INTERPOLATION REL- C ATIVE TO ARC-LENGTH ALONG THE GEODESIC. THE TANGENTIAL C COMPONENT OF G IS THE DERIVATIVE (WITH RESPECT TO ARC- C LENGTH) OF THE CUBIC INTERPOLANT AT P, WHILE THE NORMAL C COMPONENT OF G IS OBTAINED BY LINEAR INTERPOLATION OF THE C NORMAL COMPONENTS OF THE GRADIENTS AT P1 AND P2. THIS C ALGORITHM IS DUE TO C. L. LAWSON. C C INPUT PARAMETERS - P - CARTESIAN COORDINATES OF A POINT C LYING ON THE ARC DEFINED BY P1 AND C P2. C C P1,P2 - COORDINATES OF DISTINCT POINTS ON C THE UNIT SPHERE DEFINING AN ARC C WITH LENGTH LESS THAN 180 DEGREES. C C W1,W2 - DATA VALUES ASSOCIATED WITH P1 AND C P2, RESPECTIVELY. C C G1,G2 - GRADIENT VECTORS ASSOCIATED WITH P1 C AND P2. G1 AND G2 ARE ORTHOGONAL C TO P1 AND P2, RESPECTIVELY. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C G - ARRAY OF LENGTH 3. C C OUTPUT PARAMETERS - W - INTERPOLATED VALUE AT P. C C G - INTERPOLATED GRADIENT AT P. C C GN - NORMAL COMPONENT OF G WITH THE C DIRECTION P1 X P2 TAKEN TO BE C POSITIVE. THE EXTRAPOLATION C PROCEDURE REQUIRES THIS COMPONENT. C C FOR EACH VECTOR V, V(1), V(2), AND V(3) CONTAIN X-, Y-, C AND Z-COMPONENTS, RESPECTIVELY. C C MODULES REFERENCED BY ARCINT - ARCLEN C C INTRINSIC FUNCTION CALLED BY ARCINT - SQRT C C*********************************************************** C INTEGER I, LUN REAL UN(3), UNORM, TAU1, TAU2, A, AL, S, T, GT DATA LUN/6/ C C LOCAL PARAMETERS - C C I = DO-LOOP INDEX C LUN = LOGICAL UNIT FOR ERROR MESSAGES C UN = UNIT NORMAL TO THE PLANE OF P, P1, AND P2 C UNORM = EUCLIDEAN NORM OF P1 X P2 -- USED TO NORMALIZE C UN C TAU1,TAU2 = TANGENTIAL DERIVATIVES (COMPONENTS OF G1,G2) C AT P1 AND P2 C A = ANGLE IN RADIANS (ARC-LENGTH) BETWEEN P1 AND C P2 C AL = ARC-LENGTH BETWEEN P1 AND P C S = NORMALIZED VALUE OF AL -- AS P VARIES FROM P1 C TO P2, S VARIES FROM 0 TO 1 C T = 1-S -- S AND T ARE BARYCENTRIC COORDINATES OF C P WITH RESPECT TO THE ARC FROM P1 TO P2 C GT = TANGENTIAL COMPONENT OF G -- COMPONENT IN THE C DIRECTION UN X P C C C COMPUTE UNIT NORMAL UN C UN(1) = P1(2)*P2(3) - P1(3)*P2(2) UN(2) = P1(3)*P2(1) - P1(1)*P2(3) UN(3) = P1(1)*P2(2) - P1(2)*P2(1) UNORM = SQRT(UN(1)*UN(1) + UN(2)*UN(2) + UN(3)*UN(3)) IF (UNORM .EQ. 0.) GO TO 2 C C NORMALIZE UN C DO 1 I = 1,3 1 UN(I) = UN(I)/UNORM C C COMPUTE TANGENTIAL DERIVATIVES AT THE ENDPOINTS -- C TAU1 = (G1,UN X P1) = (G1,P2)/UNORM AND C TAU2 = (G2,UN X P2) = -(G2,P1)/UNORM. C TAU1 = (G1(1)*P2(1) + G1(2)*P2(2) + G1(3)*P2(3))/UNORM TAU2 =-(G2(1)*P1(1) + G2(2)*P1(2) + G2(3)*P1(3))/UNORM C C COMPUTE ARC-LENGTHS A, AL C A = ARCLEN(P1,P2) IF (A .EQ. 0.) GO TO 2 AL = ARCLEN(P1,P) C C COMPUTE W BY HERMITE CUBIC INTERPOLATION C S = AL/A T = 1. - S W = W1*(2.*S+1.)*T*T + W2*(3.-2.*S)*S*S + . A*S*T*(TAU1*T - TAU2*S) C C COMPUTE TANGENTIAL AND NORMAL DERIVATIVES AT P C GT = 6.*S*T*(W2-W1)/A + . TAU1*T*(1.-3.*S) + TAU2*S*(3.*S-2.) GN = T*(UN(1)*G1(1) + UN(2)*G1(2) + UN(3)*G1(3)) + . S*(UN(1)*G2(1) + UN(2)*G2(2) + UN(3)*G2(3)) C C COMPUTE G = GT*(UN X P) + GN*UN C G(1) = GT*(UN(2)*P(3) - UN(3)*P(2)) + GN*UN(1) G(2) = GT*(UN(3)*P(1) - UN(1)*P(3)) + GN*UN(2) G(3) = GT*(UN(1)*P(2) - UN(2)*P(1)) + GN*UN(3) RETURN C C P1 X P2 = 0. PRINT AN ERROR MESSAGE AND TERMINATE C PROCESSING. C 2 WRITE (LUN,100) (P1(I),I=1,3), (P2(I),I=1,3) 100 FORMAT (1H1,24HERROR IN ARCINT -- P1 = ,2(F9.6,3H, ), . F9.6/1H ,19X,5HP2 = ,2(F9.6,3H, ),F9.6) STOP END REAL FUNCTION ARCLEN (P,Q) REAL P(3), Q(3) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS FUNCTION COMPUTES THE ARC-LENGTH (ANGLE IN RADIANS) C BETWEEN A PAIR OF POINTS ON THE UNIT SPHERE. C C INPUT PARAMETERS - P,Q - VECTORS OF LENGTH 3 CONTAINING C THE X-, Y-, AND Z-COORDINATES (IN C THAT ORDER) OF POINTS ON THE UNIT C SPHERE. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - ARCLEN - ANGLE IN RADIANS BETWEEN THE C UNIT VECTORS P AND Q. 0 .LE. C ARCLEN .LE. PI. C C MODULES REFERENCED BY ARCLEN - NONE C C INTRINSIC FUNCTIONS CALLED BY ARCLEN - ATAN, SQRT C C*********************************************************** C INTEGER I REAL D C C LOCAL PARAMETERS - C C I = DO-LOOP INDEX C D = EUCLIDEAN NORM SQUARED OF P+Q C D = 0. DO 1 I = 1,3 1 D = D + (P(I) + Q(I))**2 IF (D .EQ. 0.) GO TO 2 IF (D .GE. 4.) GO TO 3 ARCLEN = 2.*ATAN(SQRT((4.-D)/D)) RETURN C C P AND Q ARE SEPARATED BY 180 DEGREES C 2 ARCLEN = 4.*ATAN(1.) RETURN C C P AND Q COINCIDE C 3 ARCLEN = 0. RETURN END SUBROUTINE BDYADD (KK,I1,I2, IADJ,IEND ) INTEGER KK, I1, I2, IADJ(1), IEND(KK) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE ADDS A BOUNDARY NODE TO A TRIANGULATION OF C A SET OF KK-1 POINTS ON THE UNIT SPHERE. IADJ AND IEND C ARE UPDATED WITH THE INSERTION OF NODE KK. C C INPUT PARAMETERS - KK - INDEX OF AN EXTERIOR NODE TO BE C ADDED. KK .GE. 4. C C I1 - FIRST (RIGHTMOST AS VIEWED FROM C KK) BOUNDARY NODE IN THE MESH C WHICH IS VISIBLE FROM KK - THE C LINE SEGMENT KK-I1 INTERSECTS C NO ARCS. C C I2 - LAST (LEFTMOST) BOUNDARY NODE C WHICH IS VISIBLE FROM KK. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IADJ AND IEND MAY BE CREATED BY TRMESH AND MUST CONTAIN C THE VERTICES I1 AND I2. I1 AND I2 MAY BE DETERMINED BY C TRFIND. C C KK, I1, AND I2 ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION C OF NODE KK AS THE LAST C ENTRY. NODE KK WILL BE C CONNECTED TO I1, I2, AND C ALL BOUNDARY NODES BETWEEN C THEM. NO OPTIMIZATION OF C THE MESH IS PERFORMED. C C MODULE REFERENCED BY BDYADD - SHIFTD C C INTRINSIC FUNCTIONS CALLED BY BDYADD - MIN0, MAX0 C C*********************************************************** C INTEGER K, KM1, NRIGHT, NLEFT, NF, NL, N1, N2, I, . IMIN, IMAX, KEND, NEXT, INDX C C LOCAL PARAMETERS - C C K = LOCAL COPY OF KK C KM1 = K - 1 C NRIGHT,NLEFT = LOCAL COPIES OF I1, I2 C NF,NL = INDICES OF IADJ BOUNDING THE PORTION OF THE C ARRAY TO BE SHIFTED C N1 = IADJ INDEX OF THE FIRST NEIGHBOR OF NLEFT C N2 = IADJ INDEX OF THE LAST NEIGHBOR OF NRIGHT C I = DO-LOOP INDEX C IMIN,IMAX = BOUNDS ON DO-LOOP INDEX -- FIRST AND LAST C ELEMENTS OF IEND TO BE INCREMENTED C KEND = POINTER TO THE LAST NEIGHBOR OF K IN IADJ C NEXT = NEXT BOUNDARY NODE TO BE CONNECTED TO KK C INDX = INDEX FOR IADJ C K = KK KM1 = K - 1 NRIGHT = I1 NLEFT = I2 C C INITIALIZE VARIABLES C NL = IEND(KM1) N1 = 1 IF (NLEFT .NE. 1) N1 = IEND(NLEFT-1) + 1 N2 = IEND(NRIGHT) NF = MAX0(N1,N2) C C INSERT K AS A NEIGHBOR OF MAX(NRIGHT,NLEFT) C CALL SHIFTD(NF,NL,2, IADJ) IADJ(NF+1) = K IMIN = MAX0(NRIGHT,NLEFT) DO 1 I = IMIN,KM1 1 IEND(I) = IEND(I) + 2 C C INITIALIZE KEND AND INSERT K AS A NEIGHBOR OF C MIN(NRIGHT,NLEFT) C KEND = NL + 3 NL = NF - 1 NF = MIN0(N1,N2) CALL SHIFTD(NF,NL,1, IADJ) IADJ(NF) = K IMAX = IMIN - 1 IMIN = MIN0(NRIGHT,NLEFT) DO 2 I = IMIN,IMAX 2 IEND(I) = IEND(I) + 1 C C INSERT NRIGHT AS THE FIRST NEIGHBOR OF K C IADJ(KEND) = NRIGHT C C INITIALIZE INDX FOR LOOP ON BOUNDARY NODES BETWEEN NRIGHT C AND NLEFT C INDX = IEND(NRIGHT) - 2 3 NEXT = IADJ(INDX) IF (NEXT .EQ. NLEFT) GO TO 4 C C CONNECT NEXT AND K C KEND = KEND + 1 IADJ(KEND) = NEXT INDX = IEND(NEXT) IADJ(INDX) = K INDX = INDX - 1 GO TO 3 C C INSERT NLEFT AND 0 AS THE LAST NEIGHBORS OF K C 4 IADJ(KEND+1) = NLEFT KEND = KEND + 2 IADJ(KEND) = 0 IEND(K) = KEND RETURN END SUBROUTINE BNODES (N,IADJ,IEND, NB,NA,NT,NODES) INTEGER N, IADJ(1), IEND(N), NB, NA, NT, NODES(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF N POINTS ON THE UNIT SPHERE, C THIS SUBROUTINE RETURNS A VECTOR CONTAINING THE INDICES C (IF ANY) OF THE COUNTERCLOCKWISE-ORDERED SEQUENCE OF NODES C ON THE BOUNDARY OF THE CONVEX HULL OF THE SET OF POINTS. C THE BOUNDARY IS EMPTY IF THE POINTS DO NOT LIE IN A SINGLE C HEMISPHERE. THE NUMBERS OF BOUNDARY NODES, ARCS, AND C TRIANGLES ARE ALSO RETURNED. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C C IADJ - SET OF ADJACENCY LISTS OF C NODES IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C NODES - VECTOR OF LENGTH .GE. NB. C (NB .LE. N). C C IADJ AND IEND MAY BE CREATED BY TRMESH AND ARE NOT C ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - NB - NUMBER OF BOUNDARY NODES. C C NA,NT - NUMBER OF ARCS AND TRIANGLES, C RESPECTIVELY, IN THE MESH. C C NODES - VECTOR OF NB BOUNDARY NODE C INDICES RANGING FROM 1 TO N. C C MODULES REFERENCED BY BNODES - NONE C C*********************************************************** C INTEGER NN, NST, INDL, K, N0, INDF C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C NST = FIRST ELEMENT OF NODES -- ARBITRARILY CHOSEN C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF NST C K = NODES INDEX C N0 = BOUNDARY NODE TO BE ADDED TO NODES C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF N0 C NN = N C C SEARCH FOR A BOUNDARY NODE C DO 1 NST = 1,NN INDL = IEND(NST) IF (IADJ(INDL) .EQ. 0) GO TO 2 1 CONTINUE C C NO BOUNDARY NODE EXISTS C NB = 0 NT = 2*(NN-2) NA = 3*(NN-2) RETURN C C NST IS THE FIRST BOUNDARY NODE ENCOUNTERED. INITIALIZE C FOR BOUNDARY TRAVERSAL. C 2 NODES(1) = NST K = 1 N0 = NST C C TRAVERSE THE BOUNDARY IN COUNTERCLOCKWISE ORDER C 3 INDF = 1 IF (N0 .GT. 1) INDF = IEND(N0-1) + 1 N0 = IADJ(INDF) IF (N0 .EQ. NST) GO TO 4 K = K + 1 NODES(K) = N0 GO TO 3 C C TERMINATION C 4 NB = K NT = 2*(NN-1) - NB NA = 3*(NN-1) - NB RETURN END SUBROUTINE CIRCLE (N, X,Y) INTEGER N REAL X(N), Y(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE COMPUTES THE COORDINATES OF A SEQUENCE C OF EQUISPACED POINTS ON A UNIT CIRCLE. AN (N-1)-SIDED C POLYGONAL APPROXIMATION TO THE CIRCLE MAY BE PLOTTED BY C CONNECTING (X(I),Y(I)) TO (X(I+1),Y(I+1)) FOR I = 1,..., C N-1. C C INPUT PARAMETERS - N - NUMBER OF POINTS. N .GE. 2. C C X,Y - VECTORS OF LENGTH .GE. N. C C N IS NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - X,Y - COORDINATES OF POINTS ON THE C CIRCLE WHERE (X(N),Y(N)) = C (X(1),Y(1)) = (1,0). C C MODULES REFERENCED BY CIRCLE - NONE C C INTRINSIC FUNCTIONS CALLED BY CIRCLE - ATAN, FLOAT, COS, C SIN C C*********************************************************** C INTEGER NM1, I REAL DTH, TH C C LOCAL PARAMETERS - C C NM1 = N - 1 C I = DO-LOOP INDEX C DTH = ANGLE BETWEEN ADJACENT POINTS C TH = POLAR COORDINATE ANGLE C NM1 = N - 1 IF (NM1 .LT. 1) RETURN DTH = 8.*ATAN(1.)/FLOAT(NM1) DO 1 I = 1,NM1 TH = FLOAT(I-1)*DTH X(I) = COS(TH) 1 Y(I) = SIN(TH) X(N) = X(1) Y(N) = Y(1) RETURN END SUBROUTINE CONSTR (XK,YK,ZK, CX,SX,CY,SY) REAL XK, YK, ZK, CX, SX, CY, SY C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE CONSTRUCTS THE ELEMENTS OF A 3 BY 3 C ORTHOGONAL MATRIX R WHICH ROTATES A POINT (XK,YK,ZK) ON C THE UNIT SPHERE TO THE NORTH POLE, I.E. C C (XK) (CY 0 -SY) (1 0 0) (XK) (0) C R * (YK) = ( 0 1 0) * (0 CX -SX) * (YK) = (0) C (ZK) (SY 0 CY) (0 SX CX) (ZK) (1) C C INPUT PARAMETERS - XK,YK,ZK - COMPONENTS OF A UNIT VECTOR C TO BE ROTATED TO (0,0,1). C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - CX,SX,CY,SY - ELEMENTS OF R -- CX,SX C DEFINE A ROTATION ABOUT C THE X-AXIS AND CY,SY DE- C FINE A ROTATION ABOUT C THE Y-AXIS. C C MODULES REFERENCED BY CONSTR - NONE C C INTRINSIC FUNCTION CALLED BY CONSTR - SQRT C C*********************************************************** C CY = SQRT(YK*YK + ZK*ZK) SY = XK IF (CY .EQ. 0.) GO TO 1 CX = ZK/CY SX = YK/CY RETURN C C (XK,YK,ZK) LIES ON THE X-AXIS C 1 CX = 1. SX = 0. RETURN END SUBROUTINE COVSPH (KK,NODE, IADJ,IEND ) INTEGER KK, NODE, IADJ(1), IEND(KK) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE CONNECTS AN EXTERIOR NODE KK TO ALL C BOUNDARY NODES OF A TRIANGULATION OF KK-1 POINTS ON THE C UNIT SPHERE, PRODUCING A TRIANGULATION WHICH COVERS THE C SPHERE. IADJ AND IEND ARE UPDATED WITH THE ADDITION OF C NODE KK, BUT NO OPTIMIZATION OF THE MESH IS PERFORMED. C ALL BOUNDARY NODES MUST BE VISIBLE FROM KK. C C INPUT PARAMETERS - KK - INDEX OF THE EXTERIOR NODE TO C BE ADDED. KK .GE. 4. C C NODE - BOUNDARY NODE INDEX IN THE C RANGE 1,...,KK-1. C C IADJ - SET OF ADJACENCY LISTS FOR C NODES 1,...,KK-1. C C IEND - POINTERS TO THE ENDS OF ADJA- C CENCY LISTS IN IADJ FOR NODES C 1,...,KK-1. C C KK AND NODE ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION C OF NODE KK AS THE LAST C ENTRY. ALL NODES ARE C INTERIOR. C C MODULES REFERENCED BY COVSPH - NONE C C*********************************************************** C INTEGER K, ND, NEXT, KEND, INDX C C LOCAL PARAMETERS - C C K,ND = LOCAL COPIES OF KK AND NODE C NEXT = BOUNDARY NODE TO BE CONNECTED TO K C KEND = IADJ INDEX OF THE LAST NEIGHBOR OF K C INDX = IADJ INDEX C K = KK ND = NODE C C INITIALIZATION C NEXT = ND KEND = IEND(K-1) C C WALK ALONG THE BOUNDARY CONNECTING NODE K AND NEXT. K C K REPLACES 0 AS THE LAST NEIGHBOR OF NEXT. C 1 KEND = KEND + 1 IADJ(KEND) = NEXT INDX = IEND(NEXT) IADJ(INDX) = K NEXT = IADJ(INDX-1) IF (NEXT .NE. ND) GO TO 1 IEND(K) = KEND RETURN END SUBROUTINE DELETE (N,NOUT1,NOUT2, IADJ,IEND, IER) INTEGER N, NOUT1, NOUT2, IADJ(1), IEND(N), IER EXTERNAL INDEX C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE DELETES A BOUNDARY EDGE FROM A TRIANGU- C LATION OF A SET OF POINTS ON THE UNIT SPHERE. IT MAY BE C NECESSARY TO FORCE CERTAIN EDGES TO BE PRESENT BEFORE C CALLING DELETE (SEE SUBROUTINE EDGE). NOTE THAT SUBROU- C TINES EDGE, TRFIND, AND THE ROUTINES WHICH CALL TRFIND C (ADNODE, UNIF, INTRC1, AND INTRC0) SHOULD NOT BE CALLED C FOLLOWING A DELETION. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIAN- C GULATION. C C NOUT1,NOUT2 - PAIR OF ADJACENT NODES ON THE C BOUNDARY DEFINING THE ARC TO C BE REMOVED. NOUT2 MUST BE THE C LAST NONZERO NEIGHBOR OF NOUT1. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C IADJ,IEND - DATA STRUCTURE DEFINING THE C TRIANGULATION (SEE SUBROUTINE C TRMESH). C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE REMOVAL C OF THE ARC NOUT1-NOUT2 C IF IER .EQ. 0. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF NOUT1 OR NOUT2 C IS NOT ON THE C BOUNDARY. C IER = 2 IF NOUT1 OR NOUT2 C HAS ONLY 2 NONZERO C NEIGHBORS. C IER = 3 IF NOUT2 IS NOT C THE LAST NEIGHBOR C OF NOUT1. C IER = 4 IF A DELETION C WOULD DIVIDE THE C MESH INTO TWO C REGIONS. C C MODULES REFERENCED BY DELETE - SHIFTD, INDEX C C*********************************************************** C INTEGER NN, IOUT1, IOUT2, IO1, IO2, IND12, IND21, . ITEMP, IND1F, IND1L, IND2F, IND2L, NEWBD, . INDNF, INDNL, INDN0, INDFP2, INDLM3, NF, NL, . I, IMAX C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C IOUT1,IOUT2 = LOCAL COPIES OF NOUT1 AND NOUT2 C IO1,IO2 = NOUT1,NOUT2 IN ORDER OF INCREASING MAGNITUDE C IND12 = INDEX OF IO2 IN THE ADJACENCY LIST FOR IO1 C IND21 = INDEX OF IO1 IN THE ADJACENCY LIST FOR IO2 C ITEMP = TEMPORARY STORAGE LOCATION FOR PERMUTATIONS C IND1F = IADJ INDEX OF THE FIRST NEIGHBOR OF IO1 C IND1L = IADJ INDEX OF THE LAST NEIGHBOR OF IO1 C IND2F = IADJ INDEX OF THE FIRST NEIGHBOR OF IO2 C IND2L = IADJ INDEX OF THE LAST NEIGHBOR OF IO2 C NEWBD = THE NEIGHBOR COMMON TO NOUT1 AND NOUT2 C INDNF = IADJ INDEX OF THE FIRST NEIGHBOR OF NEWBD C INDNL = IADJ INDEX OF THE LAST NEIGHBOR OF NEWBD C INDN0 = INDEX OF 0 IN THE ADJACENCY LIST FOR NEWBD C BEFORE PERMUTING THE NEIGHBORS C INDFP2 = INDNF + 2 C INDLM3 = INDNL - 3 C NF,NL = BOUNDS ON THE PORTION OF IADJ TO BE SHIFTED C I = DO-LOOP INDEX C IMAX = UPPER BOUND ON DO-LOOP FOR SHIFTING IEND C NN = N IOUT1 = NOUT1 IOUT2 = NOUT2 C C INITIALIZE INDICES C IND1F = 1 IF (IOUT1 .GT. 1) IND1F = IEND(IOUT1-1) + 1 IND1L = IEND(IOUT1) IND2F = 1 IF (IOUT2 .GT. 1) IND2F = IEND(IOUT2-1) + 1 IND2L = IEND(IOUT2) NEWBD = IADJ(IND1L-2) INDN0 = INDEX(NEWBD,IOUT2,IADJ,IEND) INDNL = IEND(NEWBD) C C ORDER VERTICES SUCH THAT THE NEIGHBORS OF IO1 PRECEDE C THOSE OF IO2 C IF (IOUT1 .GT. IOUT2) GO TO 1 IO1 = IOUT1 IO2 = IOUT2 IND12 = IND1L - 1 IND21 = IND2F GO TO 2 1 IO1 = IOUT2 IO2 = IOUT1 IND12 = IND2F IND21 = IND1L - 1 C C CHECK FOR ERRORS C 2 IF ( (IADJ(IND1L) .NE. 0) .OR. (IADJ(IND2L) .NE. 0) ) . GO TO 21 IF ( (IND1L-IND1F .LE. 2) .OR. (IND2L-IND2F .LE. 2) ) . GO TO 22 IF (IADJ(IND1L-1) .NE. IOUT2) GO TO 23 IF (IADJ(INDNL) .EQ. 0) GO TO 24 C C DELETE THE EDGE IO1-IO2 AND MAKE NEWBD A BOUNDARY NODE C IF (NEWBD .LT. IO1) GO TO 8 IF (NEWBD .LT. IO2) GO TO 6 C C THE VERTICES ARE ORDERED IO1, IO2, NEWBD. C DELETE IO2 AS A NEIGHBOR OF IO1. C NF = IND12 + 1 NL = IND21 - 1 CALL SHIFTD(NF,NL,-1, IADJ) IMAX = IO2 - 1 DO 3 I = IO1,IMAX 3 IEND(I) = IEND(I) - 1 C C DELETE IO1 AS A NEIGHBOR OF IO2 C NF = NL + 2 NL = INDN0 CALL SHIFTD(NF,NL,-2, IADJ) IMAX = NEWBD - 1 DO 4 I = IO2,IMAX 4 IEND(I) = IEND(I) - 2 C C SHIFT THE BOTTOM OF IADJ UP 1 LEAVING ROOM FOR 0 AS A C NEIGHBOR OF NEWBD C INDN0 = INDN0 - 1 NF = NL + 1 NL = IEND(NN) IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ) DO 5 I = NEWBD,NN 5 IEND(I) = IEND(I) - 1 GO TO 12 C C THE VERTICES ARE ORDERED IO1, NEWBD, IO2. C DELETE IO2 AS A NEIGHBOR OF IO1 LEAVING ROOM FOR 0 AS A C NEIGHBOR OF NEWBD. C 6 NF = IND12 + 1 NL = INDN0 CALL SHIFTD(NF,NL,-1, IADJ) IMAX = NEWBD - 1 DO 7 I = IO1,IMAX 7 IEND(I) = IEND(I) - 1 GO TO 10 C C THE VERTICES ARE ORDERED NEWBD, IO1, IO2. C DELETE IO2 AS A NEIGHBOR OF IO1 LEAVING ROOM FOR 0 AS A C NEIGHBOR OF NEWBD. C 8 INDN0 = INDN0 + 1 NF = INDN0 NL = IND12 - 1 IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ) IMAX = IO1 - 1 DO 9 I = NEWBD,IMAX 9 IEND(I) = IEND(I) + 1 C C DELETE IO1 AS A NEIGHBOR OF IO2 C 10 NF = IND21 + 1 NL = IEND(NN) CALL SHIFTD(NF,NL,-1, IADJ) DO 11 I = IO2,NN 11 IEND(I) = IEND(I) - 1 C C PERMUTE THE NEIGHBORS OF NEWBD WITH END-AROUND SHIFTS SO C THAT 0 IS THE LAST NEIGHBOR C 12 INDNF = 1 IF (NEWBD .GT. 1) INDNF = IEND(NEWBD-1) + 1 INDNL = IEND(NEWBD) IF (INDN0-INDNF .GE. INDNL-INDN0) GO TO 16 C C SHIFT UPWARD C IF (INDN0 .GT. INDNF) GO TO 13 CALL SHIFTD(INDNF+1,INDNL,-1, IADJ) GO TO 20 13 INDFP2 = INDNF + 2 IF (INDN0 .LT. INDFP2) GO TO 15 DO 14 I = INDFP2,INDN0 ITEMP = IADJ(INDNF) CALL SHIFTD(INDNF+1,INDNL,-1, IADJ) 14 IADJ(INDNL) = ITEMP C C THE LAST SHIFT IS BY 2 C 15 ITEMP = IADJ(INDNF) CALL SHIFTD(INDFP2,INDNL,-2, IADJ) IADJ(INDNL-1) = ITEMP GO TO 20 C C SHIFT DOWNWARD C 16 IF (INDN0 .EQ. INDNL) GO TO 20 IF (INDN0 .LT. INDNL-1) GO TO 17 CALL SHIFTD(INDNF,INDNL-2,1, IADJ) IADJ(INDNF) = IADJ(INDNL) GO TO 20 17 INDLM3 = INDNL - 3 IF (INDN0 .GT. INDLM3) GO TO 19 DO 18 I = INDN0,INDLM3 ITEMP = IADJ(INDNL) CALL SHIFTD(INDNF,INDNL-1,1, IADJ) 18 IADJ(INDNF) = ITEMP C C THE LAST SHIFT IS BY 2 C 19 ITEMP = IADJ(INDNL-1) CALL SHIFTD(INDNF,INDLM3,2, IADJ) IADJ(INDNF+1) = IADJ(INDNL) IADJ(INDNF) = ITEMP C C INSERT 0 AS THE LAST NEIGHBOR OF NEWBD C 20 IADJ(INDNL) = 0 IER = 0 RETURN C C ONE OF THE VERTICES IS NOT ON THE BOUNDARY C 21 IER = 1 RETURN C C ONE OF THE VERTICES HAS ONLY TWO NONZERO NEIGHBORS -- THE C TRIANGULATION WOULD BE DESTROYED BY A DELETION C 22 IER = 2 RETURN C C NOUT2 IS NOT THE LAST NONZERO NEIGHBOR OF NOUT1 C 23 IER = 3 RETURN C C A DELETION WOULD DIVIDE THE MESH INTO TWO REGIONS C CONNECTED AT A SINGLE NODE C 24 IER = 4 RETURN END SUBROUTINE EDGE (IN1,IN2,X,Y,Z, LWK,IWK,IADJ, . IEND, IER) LOGICAL SWPTST INTEGER IN1, IN2, LWK, IWK(2,1), IADJ(1), IEND(1), IER REAL X(1), Y(1), Z(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF N NODES AND A PAIR OF NODAL C INDICES IN1 AND IN2, THIS ROUTINE SWAPS ARCS AS NECESSARY C TO FORCE IN1 AND IN2 TO BE ADJACENT. ONLY ARCS WHICH C INTERSECT IN1-IN2 ARE SWAPPED OUT. IF A THIESSEN TRIANGU- C LATION IS INPUT, THE RESULTING TRIANGULATION IS AS CLOSE C AS POSSIBLE TO A THIESSEN TRIANGULATION IN THE SENSE THAT C ALL ARCS OTHER THAN IN1-IN2 ARE LOCALLY OPTIMAL. C A SEQUENCE OF CALLS TO EDGE MAY BE USED TO FORCE THE C PRESENCE OF A SET OF EDGES DEFINING THE BOUNDARY OF A NON- C CONVEX REGION. SUBSEQUENT DELETION OF EDGES OUTSIDE THIS C REGION (BY SUBROUTINE DELETE) RESULTS IN A NONCONVEX TRI- C ANGULATION WHICH MAY SERVE AS A FINITE ELEMENT GRID. C (EDGE SHOULD NOT BE CALLED AFTER A CALL TO DELETE.) IF, C ON THE OTHER HAND, INTERPOLATION IS TO BE PERFORMED IN THE C NONCONVEX REGION, EDGES MUST NOT BE DELETED, BUT IT IS C STILL ADVANTAGEOUS TO HAVE THE NONCONVEX BOUNDARY PRESENT C IF IT IS DESIRABLE THAT INTERPOLATED VALUES BE INFLUENCED C BY THE GEOMETRY. NOTE THAT SUBROUTINE GETNP WHICH IS USED C TO SELECT THE NODES ENTERING INTO LOCAL DERIVATIVE ESTI- C MATES WILL NOT NECESSARILY RETURN CLOSEST NODES IF THE C TRIANGULATION HAS BEEN RENDERED NONOPTIMAL BY A CALL TO C EDGE. HOWEVER, THE EFFECT WILL BE MERELY TO FURTHER EN- C HANCE THE INFLUENCE OF THE NONCONVEX GEOMETRY ON INTERPO- C LATED VALUES. C C INPUT PARAMETERS - IN1,IN2 - INDICES (OF X,Y,Z) IN THE C RANGE 1,...,N DEFINING A PAIR C OF NODES TO BE CONNECTED BY C AN ARC. C C X,Y,Z - N-VECTORS CONTAINING CARTE- C SIAN COORDINATES OF THE C NODES. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C LWK - NUMBER OF COLUMNS RESERVED C FOR IWK. THIS MUST BE AT C LEAST NI -- THE NUMBER OF C ARCS WHICH INTERSECT IN1-IN2. C (NI IS BOUNDED BY N-3). C C IWK - INTEGER WORK ARRAY DIMENSION- C ED 2 BY LWK (OR VECTOR OF C LENGTH .GE. 2*LWK). C C IADJ,IEND - DATA STRUCTURE DEFINING THE C TRIANGULATION. SEE SUBROU- C TINE TRMESH. C C OUTPUT PARAMETERS - LWK - NUMBER OF IWK COLUMNS REQUIRED C IF IER = 0 OR IER = 2. LWK = 0 C IFF IN1 AND IN2 WERE ADJACENT C ON INPUT. C C IWK - CONTAINS THE INDICES OF THE END- C POINTS OF THE NEW ARCS OTHER C THAN IN1-IN2 UNLESS IER .GT. 0 C OR LWK = 0. NEW ARCS TO THE C LEFT OF IN1->IN2 ARE STORED IN C THE FIRST K-1 COLUMNS (LEFT POR- C TION OF IWK), COLUMN K CONTAINS C ZEROS, AND NEW ARCS TO THE RIGHT C OF IN1->IN2 OCCUPY COLUMNS K+1, C ...,LWK. (K CAN BE DETERMINED C BY SEARCHING IWK FOR THE ZEROS.) C C IADJ,IEND - UPDATED IF NECESSARY TO REFLECT C THE PRESENCE OF AN ARC CONNECT- C ING IN1 AND IN2, UNALTERED IF C IER .NE. 0. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE EN- C COUNTERED. C IER = 1 IF IN1 .LT. 1, IN2 .LT. C 1, IN1 = IN2, OR LWK C .LT. 0 ON INPUT. C IER = 2 IF MORE SPACE IS REQUIR- C ED IN IWK. SEE LWK. C IER = 3 IF IN1 AND IN2 COULD NOT C BE CONNECTED DUE TO AN C INVALID DATA STRUCTURE. C C MODULES REFERENCED BY EDGE - SWAP, INDEX, SHIFTD, SWPTST C C*********************************************************** C INTEGER N1, N2, IWEND, IWL, INDF, INDX, N1LST, NL, NR, . NEXT, IWF, LFT, N0, IWC, IWCP1, IWCM1, I, IO1, . IO2, INDL REAL X1, Y1, Z1, X2, Y2, Z2, X0, Y0, Z0 LOGICAL SWP, LEFT C C LOCAL PARAMETERS - C C N1,N2 = LOCAL COPIES OF IN1 AND IN2 OR NODES OPPOSITE C AN ARC IO1-IO2 TO BE TESTED FOR A SWAP IN C THE OPTIMIZATION LOOPS C IWEND = INPUT OR OUTPUT VALUE OF LWK C IWL = IWK (COLUMN) INDEX OF THE LAST (RIGHTMOST) ARC C WHICH INTERSECTS IN1->IN2 C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF IN1 OR IO1 C INDX = IADJ INDEX OF A NEIGHBOR OF IN1, NL, OR IO1 C N1LST = LAST NEIGHBOR OF IN1 C NL,NR = ENDPOINTS OF AN ARC WHICH INTERSECTS IN1-IN2 C WITH NL LEFT IN1->IN2 C NEXT = NODE OPPOSITE NL->NR C IWF = IWK (COLUMN) INDEX OF THE FIRST (LEFTMOST) ARC C WHICH INTERSECTS IN1->IN2 C LFT = FLAG USED TO DETERMINE IF A SWAP RESULTS IN THE C NEW ARC INTERSECTING IN1-IN2 -- LFT = 0 IFF C N0 = IN1, LFT = -1 IMPLIES N0 LEFT IN1->IN2, C AND LFT = 1 IMPLIES N0 LEFT IN2->IN1 C N0 = NODE OPPOSITE NR->NL C IWC = IWK INDEX BETWEEN IWF AND IWL -- NL->NR IS C STORED IN IWK(1,IWC)->IWK(2,IWC) C IWCP1 = IWC + 1 C IWCM1 = IWC - 1 C I = DO-LOOP INDEX AND COLUMN INDEX FOR IWK C IO1,IO2 = ENDPOINTS OF AN ARC TO BE TESTED FOR A SWAP IN C THE OPTIMIZATION LOOPS C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF IO1 C X1,Y1,Z1 = COORDINATES OF IN1 C X2,Y2,Z2 = COORDINATES OF IN2 C X0,Y0,Z0 = COORDINATES OF N0 C SWP = FLAG SET TO .TRUE. IFF A SWAP OCCURS IN AN OPT- C IMIZATION LOOP C LEFT = STATEMENT FUNCTION WHICH RETURNS THE VALUE C .TRUE. IFF (XP,YP,ZP) IS ON OR TO THE LEFT OF C THE VECTOR (XA,YA,ZA)->(XB,YB,ZB) C LEFT(XA,YA,ZA,XB,YB,ZB,XP,YP,ZP) = XP*(YA*ZB-YB*ZA) . - YP*(XA*ZB-XB*ZA) + ZP*(XA*YB-XB*YA) .GE. 0. C C STORE IN1, IN2, AND LWK IN LOCAL VARIABLES AND CHECK FOR C ERRORS. C N1 = IN1 N2 = IN2 IWEND = LWK IF (N1 .LT. 1 .OR. N2 .LT. 1 .OR. N1 .EQ. N2 .OR. . IWEND .LT. 0) GO TO 35 C C STORE THE COORDINATES OF N1 AND N2 AND INITIALIZE IWL. C X1 = X(N1) Y1 = Y(N1) Z1 = Z(N1) X2 = X(N2) Y2 = Y(N2) Z2 = Z(N2) IWL = 0 C C SET NR AND NL TO ADJACENT NEIGHBORS OF N1 SUCH THAT C NR LEFT N2->N1 AND NL LEFT N1->N2. C C SET INDF AND INDX TO THE INDICES OF THE FIRST AND LAST C NEIGHBORS OF N1 AND SET N1LST TO THE LAST NEIGHBOR. C INDF = 1 IF (N1 .GT. 1) INDF = IEND(N1-1) + 1 INDX = IEND(N1) N1LST = IADJ(INDX) IF (N1LST .EQ. 0) INDX = INDX - 1 IF (N1LST .EQ. 0) GO TO 2 C C N1 IS AN INTERIOR NODE. LOOP THROUGH THE NEIGHBORS NL C IN REVERSE ORDER UNTIL NL LEFT N1->N2. C NL = N1LST 1 IF ( LEFT(X1,Y1,Z1,X2,Y2,Z2,X(NL),Y(NL),Z(NL)) ) . GO TO 2 INDX = INDX - 1 NL = IADJ(INDX) IF (INDX .GT. INDF) GO TO 1 C C NL IS THE FIRST NEIGHBOR OF N1. SET NR TO THE LAST C NEIGHBOR AND TEST FOR AN ARC N1-N2. C NR = N1LST IF (NL .EQ. N2) GO TO 34 GO TO 4 C C NL = IADJ(INDX) LEFT N1->N2 AND INDX .GT. INDF. SET C NR TO THE PRECEDING NEIGHBOR OF N1. C 2 INDX = INDX - 1 NR = IADJ(INDX) IF ( LEFT(X2,Y2,Z2,X1,Y1,Z1,X(NR),Y(NR),Z(NR)) ) . GO TO 3 IF (INDX .GT. INDF) GO TO 2 C C SET NL AND NR TO THE FIRST AND LAST NEIGHBORS OF N1 AND C TEST FOR AN INVALID DATA STRUCTURE (N1 CANNOT BE A C BOUNDARY NODE AND CANNOT BE ADJACENT TO N2). C NL = NR NR = N1LST IF (NR .EQ. 0 .OR. NR .EQ. N2) GO TO 37 GO TO 4 C C SET NL TO THE NEIGHBOR FOLLOWING NR AND TEST FOR AN ARC C N1-N2. C 3 NL = IADJ(INDX+1) IF (NL .EQ. N2 .OR. NR .EQ. N2) GO TO 34 C C STORE THE ORDERED SEQUENCE OF INTERSECTING EDGES NL->NR IN C IWK(1,IWL)->IWK(2,IWL). C 4 IWL = IWL + 1 IF (IWL .LE. IWEND) IWK(1,IWL) = NL IF (IWL .LE. IWEND) IWK(2,IWL) = NR C C SET NEXT TO THE NEIGHBOR OF NL WHICH FOLLOWS NR. C INDX = IEND(NL) IF (IADJ(INDX) .NE. NR) GO TO 5 C C NR IS THE LAST NEIGHBOR OF NL. SET NEXT TO THE FIRST C NEIGHBOR. C INDX = 0 IF (NL .NE. 1) INDX = IEND(NL-1) GO TO 6 C C NR IS NOT THE LAST NEIGHBOR OF NL. LOOP THROUGH THE C NEIGHBORS IN REVERSE ORDER. C 5 INDX = INDX - 1 IF (IADJ(INDX) .NE. NR) GO TO 5 C C STORE NEXT, TEST FOR AN INVALID TRIANGULATION (NL->NR C CANNOT BE A BOUNDARY EDGE), AND TEST FOR TERMINATION C OF THE LOOP. C 6 NEXT = IADJ(INDX+1) IF (NEXT .EQ. 0) GO TO 37 IF (NEXT .EQ. N2) GO TO 8 C C SET NL OR NR TO NEXT. C IF ( LEFT(X1,Y1,Z1,X2,Y2,Z2,X(NEXT),Y(NEXT),Z(NEXT)) ) . GO TO 7 NR = NEXT GO TO 4 7 NL = NEXT GO TO 4 C C IWL IS THE NUMBER OF ARCS WHICH INTERSECT N1-N2. STORE C LWK AND TEST FOR SUFFICIENT SPACE. C 8 LWK = IWL IF (IWL .GT. IWEND) GO TO 36 IWEND = IWL C C INITIALIZE FOR EDGE SWAPPING LOOP -- ALL POSSIBLE SWAPS C ARE APPLIED (EVEN IF THE NEW ARC AGAIN INTERSECTS C N1-N2), ARCS TO THE LEFT OF N1->N2 ARE STORED IN THE C LEFT PORTION OF IWK, AND ARCS TO THE RIGHT ARE STORED IN C THE RIGHT PORTION. IWF AND IWL INDEX THE FIRST AND LAST C INTERSECTING ARCS. C IER = 0 IWF = 1 C C TOP OF LOOP -- SET N0 TO N1 AND NL->NR TO THE FIRST EDGE. C IWC POINTS TO THE ARC CURRENTLY BEING PROCESSED. LFT C .LE. 0 IFF N0 LEFT N1->N2. C 9 LFT = 0 N0 = N1 X0 = X1 Y0 = Y1 Z0 = Z1 NL = IWK(1,IWF) NR = IWK(2,IWF) IWC = IWF C C SET NEXT TO THE NODE OPPOSITE NL->NR UNLESS IWC IS THE C LAST ARC. C 10 IF (IWC .EQ. IWL) GO TO 21 IWCP1 = IWC + 1 NEXT = IWK(1,IWCP1) IF (NEXT .NE. NL) GO TO 15 NEXT = IWK(2,IWCP1) C C NEXT RIGHT N1->N2 AND IWC .LT. IWL. TEST FOR A POSSIBLE C SWAP. C IF ( .NOT. LEFT(X0,Y0,Z0,X(NR),Y(NR),Z(NR),X(NEXT), . Y(NEXT),Z(NEXT)) ) GO TO 13 IF (LFT .GE. 0) GO TO 11 IF ( .NOT. LEFT(X(NL),Y(NL),Z(NL),X0,Y0,Z0,X(NEXT), . Y(NEXT),Z(NEXT)) ) GO TO 13 C C REPLACE NL->NR WITH N0->NEXT. C CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND ) IWK(1,IWC) = N0 IWK(2,IWC) = NEXT GO TO 14 C C SWAP NL-NR FOR N0-NEXT, SHIFT COLUMNS IWC+1,...,IWL TO C THE LEFT, AND STORE N0-NEXT IN THE RIGHT PORTION OF C IWK. C 11 CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND ) DO 12 I = IWCP1,IWL IWK(1,I-1) = IWK(1,I) 12 IWK(2,I-1) = IWK(2,I) IWK(1,IWL) = N0 IWK(2,IWL) = NEXT IWL = IWL - 1 NR = NEXT GO TO 10 C C A SWAP IS NOT POSSIBLE. SET N0 TO NR. C 13 N0 = NR X0 = X(N0) Y0 = Y(N0) Z0 = Z(N0) LFT = 1 C C ADVANCE TO THE NEXT ARC. C 14 NR = NEXT IWC = IWC + 1 GO TO 10 C C NEXT LEFT N1->N2, NEXT .NE. N2, AND IWC .LT. IWL. C TEST FOR A POSSIBLE SWAP. C 15 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 16 IF ( .NOT. LEFT(X0,Y0,Z0,X(NR),Y(NR),Z(NR),X(NEXT), . Y(NEXT),Z(NEXT)) ) GO TO 19 C C REPLACE NL->NR WITH NEXT->N0. C CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND ) IWK(1,IWC) = NEXT IWK(2,IWC) = N0 GO TO 20 C C SWAP NL-NR FOR N0-NEXT, SHIFT COLUMNS IWF,...,IWC-1 TO C THE RIGHT, AND STORE N0-NEXT IN THE LEFT PORTION OF C IWK. C 16 CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND ) I = IWC 17 IF (I .EQ. IWF) GO TO 18 IWK(1,I) = IWK(1,I-1) IWK(2,I) = IWK(2,I-1) I = I - 1 GO TO 17 18 IWK(1,IWF) = N0 IWK(2,IWF) = NEXT IWF = IWF + 1 GO TO 20 C C A SWAP IS NOT POSSIBLE. SET N0 TO NL. C 19 N0 = NL X0 = X(N0) Y0 = Y(N0) Z0 = Z(N0) LFT = -1 C C ADVANCE TO THE NEXT ARC. C 20 NL = NEXT IWC = IWC + 1 GO TO 10 C C N2 IS OPPOSITE NL->NR (IWC = IWL). C 21 IF (N0 .EQ. N1) GO TO 24 IF (LFT .LT. 0) GO TO 22 C C N0 RIGHT N1->N2. TEST FOR A POSSIBLE SWAP. C IF ( .NOT. LEFT(X0,Y0,Z0,X(NR),Y(NR),Z(NR),X2,Y2,Z2) ) . GO TO 9 C C SWAP NL-NR FOR N0-N2 AND STORE N0-N2 IN THE RIGHT C PORTION OF IWK. C CALL SWAP(N2,N0,NL,NR, IADJ,IEND ) IWK(1,IWL) = N0 IWK(2,IWL) = N2 IWL = IWL - 1 GO TO 9 C C N0 LEFT N1->N2. TEST FOR A POSSIBLE SWAP. C 22 IF ( .NOT. LEFT(X(NL),Y(NL),Z(NL),X0,Y0,Z0,X2,Y2,Z2) ) . GO TO 9 C C SWAP NL-NR FOR N0-N2, SHIFT COLUMNS IWF,...,IWL-1 TO THE C RIGHT, AND STORE N0-N2 IN THE LEFT PORTION OF IWK. C CALL SWAP(N2,N0,NL,NR, IADJ,IEND ) 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 9 C C IWF = IWC = IWL. SWAP OUT THE LAST ARC FOR N1-N2 AND C STORE ZEROS IN IWK. C 24 CALL SWAP(N2,N1,NL,NR, IADJ,IEND ) IWK(1,IWC) = 0 IWK(2,IWC) = 0 IF (IWC .EQ. 1) GO TO 29 C C OPTIMIZATION LOOPS -- OPTIMIZE THE SET OF NEW ARCS TO THE C LEFT OF IN1->IN2. THE LOOP IS REPEATED UNTIL NO SWAPS C ARE PERFORMED. C IWCM1 = IWC - 1 25 SWP = .FALSE. DO 28 I = 1,IWCM1 IO1 = IWK(1,I) IO2 = IWK(2,I) C C SET N1 TO THE NEIGHBOR OF IO1 WHICH FOLLOWS IO2 AND SET C N2 TO THE NEIGHBOR OF IO1 WHICH PRECEDES IO2. C INDF = 1 IF (IO1 .GT. 1) INDF = IEND(IO1-1) + 1 INDL = IEND(IO1) INDX = INDL IF (IADJ(INDX) .NE. IO2) GO TO 26 C C IO2 IS THE LAST NEIGHBOR OF IO1. C N1 = IADJ(INDF) N2 = IADJ(INDX-1) GO TO 27 C C IO2 IS NOT THE LAST NEIGHBOR OF IO1. LOOP THROUGH THE C NEIGHBORS IN REVERSE ORDER. C 26 INDX = INDX - 1 IF (IADJ(INDX) .NE. IO2) GO TO 26 N1 = IADJ(INDX+1) IF (INDX .NE. INDF) N2 = IADJ(INDX-1) IF (INDX .EQ. INDF) N2 = IADJ(INDL) C C TEST IO1-IO2 FOR A SWAP. C 27 IF ( .NOT. SWPTST(N1,N2,IO1,IO2,X,Y,Z) ) GO TO 28 SWP = .TRUE. CALL SWAP(N1,N2,IO1,IO2, IADJ,IEND ) IWK(1,I) = N1 IWK(2,I) = N2 28 CONTINUE IF (SWP) GO TO 25 C C TEST FOR TERMINATION. C 29 IF (IWC .EQ. IWEND) RETURN IWCP1 = IWC + 1 C C OPTIMIZE THE SET OF NEW ARCS TO THE RIGHT OF IN1->IN2. C 30 SWP = .FALSE. DO 33 I = IWCP1,IWEND 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. C INDF = 1 IF (IO1 .GT. 1) INDF = IEND(IO1-1) + 1 INDL = IEND(IO1) INDX = INDL IF (IADJ(INDX) .NE. IO2) GO TO 31 C N1 = IADJ(INDF) N2 = IADJ(INDX-1) GO TO 32 C 31 INDX = INDX - 1 IF (IADJ(INDX) .NE. IO2) GO TO 31 N1 = IADJ(INDX+1) IF (INDX .NE. INDF) N2 = IADJ(INDX-1) IF (INDX .EQ. INDF) N2 = IADJ(INDL) C 32 IF ( .NOT. SWPTST(N1,N2,IO1,IO2,X,Y,Z) ) GO TO 33 SWP = .TRUE. CALL SWAP(N1,N2,IO1,IO2, IADJ,IEND ) IWK(1,I) = N1 IWK(2,I) = N2 33 CONTINUE IF (SWP) GO TO 30 RETURN C C IN1 AND IN2 WERE ADJACENT ON INPUT. C 34 IER = 0 LWK = 0 RETURN C C PARAMETER OUT OF RANGE C 35 IER = 1 RETURN C C INSUFFICIENT SPACE IN IWK C 36 IER = 2 RETURN C C INVALID TRIANGULATION DATA STRUCTURE C 37 IER = 3 RETURN END SUBROUTINE GETNP (X,Y,Z,IADJ,IEND,L, NPTS, DF,IER) INTEGER IADJ(1), IEND(1), L, NPTS(L), IER REAL X(1), Y(1), Z(1), DF C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A THIESSEN TRIANGULATION OF N NODES ON THE UNIT C SPHERE AND AN ARRAY NPTS CONTAINING THE INDICES OF L-1 C NODES ORDERED BY ANGULAR DISTANCE FROM NPTS(1), THIS SUB- C ROUTINE SETS NPTS(L) TO THE INDEX OF THE NEXT NODE IN THE C SEQUENCE -- THE NODE, OTHER THAN NPTS(1),...,NPTS(L-1), C WHICH IS CLOSEST TO NPTS(1). THUS, THE ORDERED SEQUENCE C OF K CLOSEST NODES TO N1 (INCLUDING N1) MAY BE DETERMINED C BY K-1 CALLS TO GETNP WITH NPTS(1) = N1 AND L = 2,3,...,K C FOR K .GE. 2. C THE ALGORITHM USES THE FACT THAT, IN A THIESSEN TRIAN- C GULATION, THE K-TH CLOSEST NODE TO A GIVEN NODE N1 IS A C NEIGHBOR OF ONE OF THE K-1 CLOSEST NODES TO N1. C C INPUT PARAMETERS - X,Y,Z - VECTORS OF LENGTH N CONTAINING C THE CARTESIAN COORDINATES OF C THE NODES. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE TRIANGULATION. C C IEND - POINTERS TO THE ENDS OF ADJA- C CENCY LISTS FOR EACH NODE IN C THE TRIANGULATION. C C L - NUMBER OF NODES IN THE SEQUENCE C ON OUTPUT. 2 .LE. L .LE. N. C C NPTS - ARRAY OF LENGTH .GE. L CONTAIN- C ING THE INDICES OF THE L-1 C CLOSEST NODES TO NPTS(1) IN THE C FIRST L-1 LOCATIONS. C C IADJ AND IEND MAY BE CREATED BY SUBROUTINE TRMESH. C C INPUT PARAMETERS OTHER THAN NPTS ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - NPTS - UPDATED WITH THE INDEX OF THE C L-TH CLOSEST NODE TO NPTS(1) IN C POSITION L UNLESS IER = 1. C C DF - INCREASING FUNCTION (NEGATIVE C COSINE) OF THE ANGULAR DISTANCE C BETWEEN NPTS(1) AND NPTS(L) C UNLESS IER = 1. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE EN- C COUNTERED. C IER = 1 IF L IS OUT OF RANGE. C C MODULES REFERENCED BY GETNP - NONE C C INTRINSIC FUNCTION CALLED BY GETNP - IABS C C*********************************************************** C INTEGER LM1, N1, I, NI, NP, INDF, INDL, INDX, NB REAL X1, Y1, Z1, DNP, DNB C C LOCAL PARAMETERS - C C LM1 = L - 1 C N1 = NPTS(1) C I = NPTS INDEX AND DO-LOOP INDEX C NI = NPTS(I) C NP = CANDIDATE FOR NPTS(L) C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF NI C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF NI C INDX = IADJ INDEX IN THE RANGE INDF,...,INDL C NB = NEIGHBOR OF NI AND CANDIDATE FOR NP C X1,Y1,Z1 = COORDINATES OF N1 C DNP,DNB = NEGATIVE COSINES OF THE ANGULAR DISTANCES FROM C N1 TO NP AND TO NB, RESPECTIVELY C LM1 = L - 1 IF (LM1 .LT. 1) GO TO 4 IER = 0 N1 = NPTS(1) X1 = X(N1) Y1 = Y(N1) Z1 = Z(N1) C C MARK THE ELEMENTS OF NPTS C DO 1 I = 1,LM1 NI = NPTS(I) 1 IEND(NI) = -IEND(NI) C C CANDIDATES FOR NP = NPTS(L) ARE THE UNMARKED NEIGHBORS C OF NODES IN NPTS. DNP IS INITIALIZED TO -COS(PI) -- C THE MAXIMUM DISTANCE. C DNP = 1. C C LOOP ON NODES NI IN NPTS C DO 2 I = 1,LM1 NI = NPTS(I) INDF = 1 IF (NI .GT. 1) INDF = IABS(IEND(NI-1)) + 1 INDL = -IEND(NI) C C LOOP ON NEIGHBORS NB OF NI C DO 2 INDX = INDF,INDL NB = IADJ(INDX) IF (NB .EQ. 0 .OR. IEND(NB) .LT. 0) GO TO 2 C C NB IS AN UNMARKED NEIGHBOR OF NI. REPLACE NP IF NB IS C CLOSER TO N1. C DNB = -(X(NB)*X1 + Y(NB)*Y1 + Z(NB)*Z1) IF (DNB .GE. DNP) GO TO 2 NP = NB DNP = DNB 2 CONTINUE NPTS(L) = NP DF = DNP C C UNMARK THE ELEMENTS OF NPTS C DO 3 I = 1,LM1 NI = NPTS(I) 3 IEND(NI) = -IEND(NI) RETURN C C L IS OUT OF RANGE C 4 IER = 1 RETURN END SUBROUTINE GIVENS ( A,B, C,S) REAL A, B, C, S C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE CONSTRUCTS THE GIVENS PLANE ROTATION -- C ( C S) C G = ( ) WHERE C*C + S*S = 1 -- WHICH ZEROS THE SECOND C (-S C) C ENTRY OF THE 2-VECTOR (A B)-TRANSPOSE. A CALL TO GIVENS C IS NORMALLY FOLLOWED BY A CALL TO ROTATE WHICH APPLIES C THE TRANSFORMATION TO A 2 BY N MATRIX. THIS ROUTINE WAS C TAKEN FROM LINPACK. C C INPUT PARAMETERS - A,B - COMPONENTS OF THE 2-VECTOR TO BE C ROTATED. C C OUTPUT PARAMETERS - A - OVERWRITTEN BY R = +/-SQRT(A*A C + B*B) C C B - OVERWRITTEN BY A VALUE Z WHICH C ALLOWS C AND S TO BE RECOVERED C AS FOLLOWS - C C = SQRT(1-Z*Z), S=Z IF ABS(Z) C .LE. 1. C C = 1/Z, S = SQRT(1-C*C) IF C ABS(Z) .GT. 1. C C C - +/-(A/R) C C S - +/-(B/R) C C MODULES REFERENCED BY GIVENS - NONE C C INTRINSIC FUNCTIONS CALLED BY GIVENS - ABS, SQRT C C*********************************************************** C REAL AA, BB, R, U, V C C LOCAL PARAMETERS - C C AA,BB = LOCAL COPIES OF A AND B C R = C*A + S*B = +/-SQRT(A*A+B*B) C U,V = VARIABLES USED TO SCALE A AND B FOR COMPUTING R C AA = A BB = B IF (ABS(AA) .LE. ABS(BB)) GO TO 1 C C ABS(A) .GT. ABS(B) C U = AA + AA V = BB/U R = SQRT(.25 + V*V) * U C = AA/R S = V * (C + C) C C NOTE THAT R HAS THE SIGN OF A, C .GT. 0, AND S HAS C SIGN(A)*SIGN(B) C B = S A = R RETURN C C ABS(A) .LE. ABS(B) C 1 IF (BB .EQ. 0.) GO TO 2 U = BB + BB V = AA/U C C STORE R IN A C A = SQRT(.25 + V*V) * U S = BB/A C = V * (S + S) C C NOTE THAT R HAS THE SIGN OF B, S .GT. 0, AND C HAS C SIGN(A)*SIGN(B) C B = 1. IF (C .NE. 0.) B = 1./C RETURN C C A = B = 0. C 2 C = 1. S = 0. RETURN END SUBROUTINE GRADG (N,X,Y,Z,W,IADJ,IEND,EPS, NIT, . GRAD, IER) INTEGER N, IADJ(1), IEND(N), NIT, IER REAL X(N), Y(N), Z(N), W(N), EPS, GRAD(3,N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF N NODES ON THE UNIT SPHERE WITH C ASSOCIATED DATA VALUES, THIS ROUTINE USES A GLOBAL METHOD C TO COMPUTE ESTIMATED GRADIENTS AT THE NODES. THE METHOD C CONSISTS OF MINIMIZING A QUADRATIC FUNCTIONAL Q(G) OVER C THE N-VECTOR G OF GRADIENTS WHERE Q APPROXIMATES THE LIN- C EARIZED CURVATURE OF AN INTERPOLANT F OVER THE TRIANGULA- C TION. THE RESTRICTION OF F TO AN ARC OF THE TRIANGULATION C IS TAKEN TO BE THE HERMITE CUBIC (WITH RESPECT TO ARC- C LENGTH) INTERPOLANT OF THE DATA VALUES AND TANGENTIAL C GRADIENT COMPONENTS AT THE ENDPOINTS OF THE ARC, AND Q IS C THE SUM OF THE LINEARIZED CURVATURES OF F ALONG THE ARCS C -- THE INTEGRALS OVER THE ARCS OF D2F(A)**2 WHERE D2F(A) C IS THE SECOND DERIVATIVE OF F WITH RESPECT TO ARC-LENGTH C A. C SINCE THE GRADIENT AT NODE K LIES IN THE PLANE TANGENT C TO THE SPHERE SURFACE AT K, IT IS DEFINED BY TWO COMPO- C NENTS -- ITS X AND Y COMPONENTS IN THE COORDINATE SYSTEM C OBTAINED BY ROTATING K TO THE NORTH POLE. THUS, THE MIN- C IMIZATION PROBLEM CORRESPONDS TO AN ORDER 2N SYMMETRIC C POSITIVE-DEFINITE SPARSE LINEAR SYSTEM WHICH IS SOLVED BY C THE BLOCK GAUSS-SEIDEL METHOD WITH 2 BY 2 BLOCKS. C AN ALTERNATIVE METHOD, SUBROUTINE GRADL, COMPUTES A C LOCAL APPROXIMATION TO THE GRADIENT AT A SINGLE NODE AND, C WHILE SLIGHTLY LESS EFFICIENT, WAS FOUND TO BE GENERALLY C MORE ACCURATE WHEN THE NODAL DISTRIBUTION IS VERY DENSE, C VARIES GREATLY, OR DOES NOT COVER THE SPHERE. GRADG, ON C THE OTHER HAND, WAS FOUND TO BE SLIGHTLY MORE ACCURATE ON C A SOMEWHAT UNIFORM DISTRIBUTION OF 514 NODES. C C INPUT PARAMETERS - N - NUMBER OF NODES. N .GE. 3. C C X,Y,Z - CARTESIAN COORDINATES OF THE NODES. C X(I)**2 + Y(I)**2 + Z(I)**2 = 1 FOR C I = 1,...,N. C C W - DATA VALUES AT THE NODES. W(I) IS C ASSOCIATED WITH (X(I),Y(I),Z(I)). C C IADJ,IEND - DATA STRUCTURE DEFINING THE TRIAN- C GULATION. SEE SUBROUTINE TRMESH. C C EPS - NONNEGATIVE CONVERGENCE CRITERION. C THE METHOD IS TERMINATED WHEN THE C MAXIMUM CHANGE IN A GRADIENT COMPO- C NENT BETWEEN ITERATIONS IS AT MOST C EPS. EPS = 1.E-3 IS SUFFICIENT FOR C EFFECTIVE CONVERGENCE. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C NIT - MAXIMUM NUMBER OF GAUSS-SEIDEL C ITERATIONS TO BE APPLIED. THIS C MAXIMUM WILL LIKELY BE ACHIEVED IF C EPS IS SMALLER THAN THE MACHINE C PRECISION. OPTIMAL EFFICIENCY WAS C ACHIEVED IN TESTING WITH EPS = 0 C AND NIT = 5 OR 6. C C GRAD - INITIAL SOLUTION ESTIMATES (ZERO C VECTORS ARE SUFFICIENT). GRAD(I,J) C CONTAINS COMPONENT I OF THE GRADI- C ENT AT NODE J FOR I = 1,2,3 (X,Y,Z) C AND J = 1,...,N. GRAD( ,J) MUST BE C ORTHOGONAL TO NODE J -- GRAD(1,J)* C X(J) + GRAD(2,J)*Y(J) + GRAD(3,J)* C Z(J) = 0. C C OUTPUT PARAMETERS - NIT - NUMBER OF GAUSS-SEIDEL ITERA- C TIONS EMPLOYED. C C GRAD - ESTIMATED GRADIENTS. SEE THE C DESCRIPTION UNDER INPUT PARAME- C TERS. GRAD IS NOT CHANGED IF C IER = 2. C C IER - ERROR INDICATOR C IER = 0 IF THE CONVERGENCE CRI- C TERION WAS ACHIEVED. C IER = 1 IF CONVERGENCE WAS NOT C ACHIEVED WITHIN NIT C ITERATIONS. C IER = 2 IF N OR EPS IS OUT OF C RANGE OR NIT .LT. 0 ON C INPUT. C C MODULES REFERENCED BY GRADG - CONSTR, APLYRT C C INTRINSIC FUNCTIONS CALLED BY GRADG - ATAN, SQRT, AMAX1, C ABS C C*********************************************************** C INTEGER NN, MAXIT, ITER, K, INDF, INDL, INDX, NB REAL TOL, DGMAX, XK, YK, ZK, WK, G1, G2, G3, CX, . SX, CY, SY, A11, A12, A22, R1, R2, XNB, YNB, . ZNB, ALFA, XS, YS, SINAL, D, T, DG1, DG2, . DGK(3) C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C MAXIT = INPUT VALUE OF NIT C ITER = NUMBER OF ITERATIONS USED C K = DO-LOOP AND NODE INDEX C INDF,INDL = IADJ INDICES OF THE FIRST AND LAST NEIGHBORS C OF K C INDX = IADJ INDEX IN THE RANGE INDF,...,INDL C NB = NEIGHBOR OF K C TOL = LOCAL COPY OF EPS C DGMAX = MAXIMUM CHANGE IN A GRADIENT COMPONENT BE- C TWEEN ITERATIONS C XK,YK,ZK,WK = X(K), Y(K), Z(K), W(K) C G1,G2,G3 = COMPONENTS OF GRAD( ,K) C CX,SX,CY,SY = COMPONENTS OF A ROTATION MAPPING K TO THE C NORTH POLE (0,0,1) C A11,A12,A22 = MATRIX COMPONENTS OF THE 2 BY 2 BLOCK A*DG C = R WHERE A IS SYMMETRIC, (DG1,DG2,0) IS C THE CHANGE IN THE GRADIENT AT K, AND R IS C THE RESIDUAL C R1,R2 = COMPONENTS OF THE RESIDUAL -- DERIVATIVES OF C Q WITH RESPECT TO THE COMPONENTS OF THE C GRADIENT AT NODE K C XNB,YNB,ZNB = COORDINATES OF NODE NB IN THE ROTATED COOR- C DINATE SYSTEM C ALFA = ARC-LENGTH BETWEEN NODES K AND NB C XS,YS = XNB**2, YNB**2 C SINAL = SIN(ALFA) -- MAGNITUDE OF THE VECTOR CROSS C PRODUCT K X NB C D = ALFA*SINAL**2 -- FACTOR IN THE 2 BY 2 SYSTEM C T = TEMPORARY STORAGE AND FACTOR OF R1 AND R2 C DG1,DG2 = SOLUTION OF THE 2 BY 2 SYSTEM -- FIRST 2 C COMPONENTS OF DGK IN THE ROTATED COORDI- C NATE SYSTEM C DGK = CHANGE IN GRAD( ,K) FROM THE PREVIOUS ESTI- C MATE C NN = N TOL = EPS MAXIT = NIT C C ERROR CHECKS AND INITIALIZATION C IF (NN .LT. 3 .OR. TOL .LT. 0. .OR. MAXIT .LT. 0) . GO TO 5 ITER = 0 C C TOP OF ITERATION LOOP C 1 IF (ITER .EQ. MAXIT) GO TO 4 DGMAX = 0. INDL = 0 DO 3 K = 1,NN XK = X(K) YK = Y(K) ZK = Z(K) WK = W(K) G1 = GRAD(1,K) G2 = GRAD(2,K) G3 = GRAD(3,K) C C CONSTRUCT THE ROTATION MAPPING NODE K TO THE NORTH POLE C CALL CONSTR (X(K),Y(K),Z(K), CX,SX,CY,SY) C C INITIALIZE COMPONENTS OF THE 2 BY 2 SYSTEM C A11 = 0. A12 = 0. A22 = 0. R1 = 0. R2 = 0. C C LOOP ON NEIGHBORS NB OF K C INDF = INDL + 1 INDL = IEND(K) DO 2 INDX = INDF,INDL NB = IADJ(INDX) IF (NB .EQ. 0) GO TO 2 C C COMPUTE THE COORDINATES OF NB IN THE ROTATED SYSTEM C T = SX*Y(NB) + CX*Z(NB) YNB = CX*Y(NB) - SX*Z(NB) ZNB = SY*X(NB) + CY*T XNB = CY*X(NB) - SY*T C C COMPUTE ARC-LENGTH ALFA BETWEN NB AND K, SINAL = C SIN(ALFA), AND D = ALFA*SIN(ALFA)**2 C ALFA = 2.*ATAN(SQRT((1.-ZNB)/(1.+ZNB))) XS = XNB*XNB YS = YNB*YNB SINAL = SQRT(XS+YS) D = ALFA*(XS+YS) C C UPDATE THE SYSTEM COMPONENTS FOR NODE NB C A11 = A11 + XS/D A12 = A12 + XNB*YNB/D A22 = A22 + YS/D T = 1.5*(W(NB)-WK)/(ALFA*ALFA*SINAL) + . ((GRAD(1,NB)*XK + GRAD(2,NB)*YK + . GRAD(3,NB)*ZK)/2. - (G1*X(NB) + . G2*Y(NB) + G3*Z(NB)))/D R1 = R1 + T*XNB R2 = R2 + T*YNB 2 CONTINUE C C SOLVE THE 2 BY 2 SYSTEM AND UPDATE DGMAX C DG2 = (A11*R2 - A12*R1)/(A11*A22 - A12*A12) DG1 = (R1 - A12*DG2)/A11 DGMAX = AMAX1(DGMAX,ABS(DG1),ABS(DG2)) C C ROTATE (DG1,DG2,0) BACK TO THE ORIGINAL COORDINATE C SYSTEM AND UPDATE GRAD( ,K) C CALL APLYRT (DG1,DG2,CX,SX,CY,SY, DGK) GRAD(1,K) = G1 + DGK(1) GRAD(2,K) = G2 + DGK(2) 3 GRAD(3,K) = G3 + DGK(3) C C INCREMENT ITER AND TEST FOR CONVERGENCE C ITER = ITER + 1 IF (DGMAX .GT. TOL) GO TO 1 C C METHOD CONVERGED C NIT = ITER IER = 0 RETURN C C METHOD FAILED TO CONVERGE WITHIN NIT ITERATIONS C 4 IER = 1 RETURN C C PARAMETER OUT OF RANGE C 5 NIT = 0 IER = 2 RETURN END SUBROUTINE GRADL (N,K,X,Y,Z,W,IADJ,IEND, G,IER) INTEGER N, K, IADJ(1), IEND(N), IER REAL X(N), Y(N), Z(N), W(N), G(3) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF A SET OF NODES ON THE UNIT C SPHERE WITH THEIR ASSOCIATED DATA VALUES W, THIS ROUTINE C ESTIMATES A GRADIENT VECTOR AT NODE K AS FOLLOWS -- THE C COORDINATE SYSTEM IS ROTATED SO THAT K BECOMES THE NORTH C POLE, NODE K AND A SET OF NEARBY NODES ARE PROJECTED C ORTHOGONALLY ONTO THE X-Y PLANE (IN THE NEW COORDINATE C SYSTEM), A QUADRATIC IS FITTED IN A WEIGHTED LEAST-SQUARES C SENSE TO THE DATA VALUES AT THE PROJECTED NODES SUCH THAT C THE VALUE (ASSOCIATED WITH K) AT (0,0) IS INTERPOLATED, X- C AND Y-PARTIAL DERIVATIVE ESTIMATES DX AND DY ARE COMPUTED C BY DIFFERENTIATING THE QUADRATIC AT (0,0), AND THE ESTI- C MATED GRADIENT G IS OBTAINED BY ROTATING (DX,DY,0) BACK TO C THE ORIGINAL COORDINATE SYSTEM. NOTE THAT G LIES IN THE C PLANE TANGENT TO THE SPHERE AT NODE K, I.E. G IS ORTHOGO- C NAL TO THE UNIT VECTOR REPRESENTED BY NODE K. A MARQUARDT C STABILIZATION FACTOR IS USED IF NECESSARY TO ENSURE A C WELL-CONDITIONED LEAST SQUARES SYSTEM, AND A UNIQUE SOLU- C TION EXISTS UNLESS THE NODES ARE COLLINEAR. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIANGU- C LATION. N .GE. 7. C C K - NODE AT WHICH THE GRADIENT IS C SOUGHT. 1 .LE. K .LE. N. C C X,Y,Z - CARTESIAN COORDINATES OF THE C NODES. C C W - DATA VALUES AT THE NODES. W(I) C IS ASSOCIATED WITH (X(I),Y(I), C Z(I)) FOR I = 1,...,N. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE TRIANGULATION. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS FOR EACH NODE. C C IADJ AND IEND MAY BE CREATED BY SUBROUTINE TRMESH. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - G - X-, Y-, AND Z-COMPONENTS (IN C THAT ORDER) OF THE ESTIMATED C GRADIENT AT NODE K UNLESS C IER .LT. 0. C C IER - ERROR INDICATOR C IER .GE. 6 IF NO ERRORS WERE C ENCOUNTERED. IER C CONTAINS THE NUMBER C OF NODES (INCLUDING C K) USED IN THE LEAST C SQUARES FIT. C IER = -1 IF N OR K IS OUT OF C RANGE. C IER = -2 IF THE LEAST SQUARES C SYSTEM HAS NO UNIQUE C SOLUTION DUE TO DUP- C LICATE OR COLLINEAR C NODES. C C MODULES REFERENCED BY GRADL - GETNP, CONSTR, APLYR, C SETUP, GIVENS, ROTATE, C APLYRT C C INTRINSIC FUNCTIONS CALLED BY GRADL - MIN0, FLOAT, SQRT, C AMIN1, ABS C C*********************************************************** C INTEGER NN, KK, LMN, LMX, LMIN, LMAX, LM1, LNP, . NPTS(30), IERR, NP, I, J, IM1, IP1, JP1, L REAL WK, SUM, DF, RF, RTOL, AVSQ, AV, RIN, CX, SX, . CY, SY, XP, YP, ZP, WT, A(6,6), C, S, DMIN, . DTOL, SF, DX, DY DATA LMN/10/ DATA LMX/30/, RTOL/1.E-6/, DTOL/.01/, SF/1./ C C LOCAL PARAMETERS - C C NN,KK = LOCAL COPIES OF N AND K C LMN,LMX = MINIMUM AND MAXIMUM VALUES OF LNP FOR N C SUFFICIENTLY LARGE. IN MOST CASES LMN-1 C NODES ARE USED IN THE FIT. 7 .LE. LMN .LE. C LMX. C LMIN,LMAX = MIN(LMN,N), MIN(LMX,N) C LM1 = LMIN-1 C LNP = LENGTH OF NPTS OR LMAX+1 C NPTS = ARRAY CONTAINING THE INDICES OF A SEQUENCE OF C NODES ORDERED BY ANGULAR DISTANCE FROM K. C NPTS(1)=K AND THE FIRST LNP-1 ELEMENTS OF C NPTS ARE USED IN THE LEAST SQUARES FIT. C UNLESS LNP = LMAX+1, NPTS(LNP) DETERMINES R C (SEE RIN). C IERR = ERROR FLAG FOR CALLS TO GETNP (NOT CHECKED) C NP = ELEMENT OF NPTS TO BE ADDED TO THE SYSTEM C I,J = LOOP INDICES C IM1,IP1 = I-1, I+1 C JP1 = J+1 C L = NUMBER OF COLUMNS OF A**T TO WHICH A ROTATION C IS APPLIED C WK = W(K) -- DATA VALUE AT NODE K C SUM = SUM OF SQUARED EUCLIDEAN DISTANCES (IN THE C ROTATED COORDINATE SYSTEM) BETWEEN THE C ORIGIN AND THE NODES USED IN THE LEAST C SQUARES FIT C DF = NEGATIVE Z-COMPONENT (IN THE ROTATED COORDI- C NATE SYSTEM) OF AN ELEMENT NP OF NPTS -- C INCREASING FUNCTION OF THE ANGULAR DISTANCE C BETWEEN K AND NP. DF LIES IN THE INTERVAL C (-1,1). C RF = VALUE OF DF ASSOCIATED WITH NPTS(LNP) UNLESS C LNP = LMAX+1 (SEE RIN) C RTOL = TOLERANCE FOR DETERMINING LNP (AND HENCE R) -- C IF THE INCREASE IN DF BETWEEN TWO SUCCESSIVE C ELEMENTS OF NPTS IS LESS THAN RTOL, THEY ARE C TREATED AS BEING THE SAME DISTANCE FROM NODE C K AND AN ADDITIONAL NODE IS ADDED C AVSQ = AV*AV -- ACCUMULATED IN SUM C AV = ROOT-MEAN-SQUARE DISTANCE (IN THE ROTATED C COORDINATE SYSTEM) BETWEEN THE ORIGIN AND C THE NODES (OTHER THAN K) IN THE LEAST C SQUARES FIT. THE FIRST 3 COLUMNS OF A**T C ARE SCALED BY 1/AVSQ, THE NEXT 2 BY 1/AV. C RIN = INVERSE OF A RADIUS OF INFLUENCE R WHICH C ENTERS INTO WT -- R = 1+RF UNLESS ALL ELE- C MENTS OF NPTS ARE USED IN THE FIT (LNP = C LMAX+1), IN WHICH CASE R IS THE DISTANCE C FUNCTION ASSOCIATED WITH SOME POINT MORE C DISTANT FROM K THAN NPTS(LMAX) C CX,SX = COMPONENTS OF A PLANE ROTATION ABOUT THE X- C AXIS WHICH, TOGETHER WITH CY AND SY, DEFINE C A MAPPING FROM NODE K TO THE NORTH POLE C (0,0,1) C CY,SY = COMPONENTS OF A PLANE ROTATION ABOUT THE Y- C AXIS C XP,YP,ZP = COORDINATES OF NP IN THE ROTATED COORDINATE C SYSTEM UNLESS ZP .LT. 0, IN WHICH CASE C (XP,YP,0) LIES ON THE EQUATOR C WT = WEIGHT FOR THE EQUATION CORRESPONDING TO NP -- C WT = (R-D)/(R*D) = 1/D - RIN WHERE D = 1-ZP C IS ASSOCIATED WITH NP C A = TRANSPOSE OF THE (UPPER TRIANGLE OF THE) AUG- C MENTED REGRESSION MATRIX C C,S = COMPONENTS OF THE PLANE ROTATION USED TO C TRIANGULARIZE THE REGRESSION MATRIX C DMIN = MINIMUM OF THE MAGNITUDES OF THE DIAGONAL C ELEMENTS OF THE TRIANGULARIZED REGRESSION C MATRIX C DTOL = TOLERANCE FOR DETECTING AN ILL-CONDITIONED C SYSTEM -- DMIN IS REQUIRED TO BE AT LEAST C DTOL C SF = MARQUARDT STABILIZATION FACTOR USED TO DAMP C OUT THE FIRST 3 SOLUTION COMPONENTS (SECOND C PARTIALS OF THE QUADRATIC) WHEN THE SYSTEM C IS ILL-CONDITIONED. INCREASING SF RESULTS C IN MORE DAMPING (A MORE NEARLY LINEAR FIT). C DX,DY = X AND Y COMPONENTS OF THE ESTIMATED GRADIENT C IN THE ROTATED COORDINATE SYSTEM C NN = N KK = K WK = W(KK) C C CHECK FOR ERRORS AND INITIALIZE LMIN, LMAX C IF (NN .LT. 7 .OR. KK .LT. 1 .OR. KK .GT. NN) . GO TO 13 LMIN = MIN0(LMN,NN) LMAX = MIN0(LMX,NN) C C COMPUTE NPTS, LNP, AVSQ, AV, AND R. C SET NPTS TO THE CLOSEST LMIN-1 NODES TO K. DF CONTAINS C THE NEGATIVE Z-COMPONENT (IN THE ROTATED COORDINATE C SYSTEM) OF THE NEW NODE ON RETURN FROM GETNP. C SUM = 0. NPTS(1) = KK LM1 = LMIN - 1 DO 1 LNP = 2,LM1 CALL GETNP (X,Y,Z,IADJ,IEND,LNP, NPTS, DF,IERR) 1 SUM = SUM + 1. - DF*DF C C ADD ADDITIONAL NODES TO NPTS UNTIL THE INCREASE IN C R = 1+RF IS AT LEAST RTOL. C DO 2 LNP = LMIN,LMAX CALL GETNP (X,Y,Z,IADJ,IEND,LNP, NPTS, RF,IERR) IF (RF-DF .GE. RTOL) GO TO 3 2 SUM = SUM + 1. - RF*RF C C USE ALL LMAX NODES IN THE LEAST SQUARES FIT. R IS C ARBITRARILY INCREASED BY 5 PER CENT. C RF = 1.05*RF + .05 LNP = LMAX + 1 C C THERE ARE LNP-2 EQUATIONS CORRESPONDING TO NODES C NPTS(2),...,NPTS(LNP-1). C 3 AVSQ = SUM/FLOAT(LNP-2) AV = SQRT(AVSQ) RIN = 1./(1.+RF) C C CONSTRUCT THE ROTATION C CALL CONSTR (X(KK),Y(KK),Z(KK), CX,SX,CY,SY) C C SET UP THE FIRST 5 EQUATIONS OF THE AUGMENTED REGRESSION C MATRIX (TRANSPOSED) AS THE COLUMNS OF A, AND ZERO OUT C THE LOWER TRIANGLE (UPPER TRIANGLE OF A) WITH GIVENS C ROTATIONS C DO 5 I = 1,5 NP = NPTS(I+1) CALL APLYR (X(NP),Y(NP),Z(NP),CX,SX,CY,SY, XP,YP,ZP) WT = 1./(1.-ZP) - RIN CALL SETUP (XP,YP,W(NP),WK,AV,AVSQ,WT, A(1,I)) IF (I .EQ. 1) GO TO 5 IM1 = I - 1 DO 4 J = 1,IM1 JP1 = J + 1 L = 6 - J CALL GIVENS ( A(J,J),A(J,I), C,S) 4 CALL ROTATE (L,C,S, A(JP1,J),A(JP1,I) ) 5 CONTINUE C C ADD THE ADDITIONAL EQUATIONS TO THE SYSTEM USING C THE LAST COLUMN OF A -- I .LE. LNP. C I = 7 6 IF (I .EQ. LNP) GO TO 8 NP = NPTS(I) CALL APLYR (X(NP),Y(NP),Z(NP),CX,SX,CY,SY, XP,YP,ZP) WT = 1./(1.-ZP) - RIN CALL SETUP (XP,YP,W(NP),WK,AV,AVSQ,WT, A(1,6)) DO 7 J = 1,5 JP1 = J + 1 L = 6 - J CALL GIVENS ( A(J,J),A(J,6), C,S) 7 CALL ROTATE (L,C,S, A(JP1,J),A(JP1,6) ) I = I + 1 GO TO 6 C C TEST THE SYSTEM FOR ILL-CONDITIONING C 8 DMIN = AMIN1( ABS(A(1,1)),ABS(A(2,2)),ABS(A(3,3)), . ABS(A(4,4)),ABS(A(5,5)) ) IF (DMIN .GE. DTOL) GO TO 12 IF (LNP .GT. LMAX) GO TO 9 C C ADD ANOTHER NODE TO THE SYSTEM AND INCREASE R -- C I .EQ. LNP C LNP = LNP + 1 IF (LNP .LE. LMAX) CALL GETNP (X,Y,Z,IADJ,IEND,LNP, . NPTS, RF,IERR) RIN = 1./(1.05*(1.+RF)) GO TO 6 C C STABILIZE THE SYSTEM BY DAMPING SECOND PARTIALS --ADD C MULTIPLES OF THE FIRST THREE UNIT VECTORS TO THE FIRST C THREE EQUATIONS. C 9 DO 11 I = 1,3 A(I,6) = SF IP1 = I + 1 DO 10 J = IP1,6 10 A(J,6) = 0. DO 11 J = I,5 JP1 = J + 1 L = 6 - J CALL GIVENS ( A(J,J),A(J,6), C,S) 11 CALL ROTATE (L,C,S, A(JP1,J),A(JP1,6) ) C C TEST THE LINEAR PORTION OF THE STABILIZED SYSTEM FOR C ILL-CONDITIONING C DMIN = AMIN1( ABS(A(4,4)),ABS(A(5,5)) ) IF (DMIN .LT. DTOL) GO TO 14 C C SOLVE THE 2 BY 2 TRIANGULAR SYSTEM FOR THE ESTIMATED C PARTIAL DERIVATIVES C 12 DY = A(6,5)/A(5,5) DX = (A(6,4) - A(5,4)*DY)/A(4,4)/AV DY = DY/AV C C ROTATE THE GRADIENT (DX,DY,0) BACK INTO THE ORIGINAL C COORDINATE SYSTEM C CALL APLYRT (DX,DY,CX,SX,CY,SY, G) IER = LNP - 1 RETURN C C N OR K IS OUT OF RANGE C 13 IER = -1 RETURN C C NO UNIQUE SOLUTION DUE TO COLLINEAR NODES C 14 IER = -2 RETURN END INTEGER FUNCTION INDEX (NVERTX,NABOR,IADJ,IEND) INTEGER NVERTX, NABOR, IADJ(1), IEND(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS FUNCTION RETURNS THE INDEX OF NABOR IN THE C ADJACENCY LIST FOR NVERTX. C C INPUT PARAMETERS - NVERTX - NODE WHOSE ADJACENCY LIST IS C TO BE SEARCHED. C C NABOR - NODE WHOSE INDEX IS TO BE C RETURNED. NABOR MUST BE C CONNECTED TO NVERTX. C C IADJ - SET OF ADJACENCY LISTS. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - INDEX - IADJ(INDEX) = NABOR. C C MODULES REFERENCED BY INDEX - NONE C C*********************************************************** C INTEGER NB, INDX C C LOCAL PARAMETERS - C C NB = LOCAL COPY OF NABOR C INDX = INDEX FOR IADJ C NB = NABOR C C INITIALIZATION C INDX = IEND(NVERTX) + 1 C C SEARCH THE LIST OF NVERTX NEIGHBORS FOR NB C 1 INDX = INDX - 1 IF (IADJ(INDX) .NE. NB) GO TO 1 C INDEX = INDX RETURN END SUBROUTINE INTADD (KK,I1,I2,I3, IADJ,IEND ) INTEGER KK, I1, I2, I3, IADJ(1), IEND(KK) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE ADDS AN INTERIOR NODE TO A TRIANGULATION C OF A SET OF KK-1 POINTS ON THE UNIT SPHERE. IADJ AND IEND C ARE UPDATED WITH THE INSERTION OF NODE KK IN THE TRIANGLE C WHOSE VERTICES ARE I1, I2, AND I3. C C INPUT PARAMETERS - KK - INDEX OF NODE TO BE C INSERTED. KK .GE. 4. C C I1,I2,I3 - INDICES OF THE VERTICES OF C A TRIANGLE CONTAINING NODE C KK - IN COUNTERCLOCKWISE C ORDER. C C IADJ - SET OF ADJACENCY LISTS C OF NODES IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IADJ AND IEND MAY BE CREATED BY TRMESH AND MUST CONTAIN C THE VERTICES I1, I2, AND I3. I1,I2,I3 MAY BE DETERMINED C BY TRFIND. C C KK, I1, I2, AND I3 ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION C OF NODE KK AS THE LAST C ENTRY. NODE KK WILL BE C CONNECTED TO NODES I1, I2, C AND I3. NO OPTIMIZATION C OF THE MESH IS PERFORMED. C C MODULE REFERENCED BY INTADD - SHIFTD C C INTRINSIC FUNCTION CALLED BY INTADD - MOD C C*********************************************************** C INTEGER K, KM1, N(3), NFT(3), IP1, IP2, IP3, INDX, NF, . NL, N1, N2, IMIN, IMAX, I, ITEMP C C LOCAL PARAMETERS - C C K = LOCAL COPY OF KK C KM1 = K - 1 C N = VECTOR CONTAINING I1, I2, I3 C NFT = POINTERS TO THE TOPS OF THE 3 SETS OF IADJ C ELEMENTS TO BE SHIFTED DOWNWARD C IP1,IP2,IP3 = PERMUTATION INDICES FOR N AND NFT C INDX = INDEX FOR IADJ AND N C NF,NL = INDICES OF FIRST AND LAST ENTRIES IN IADJ C TO BE SHIFTED DOWN C N1,N2 = FIRST 2 VERTICES OF A NEW TRIANGLE -- C (N1,N2,KK) C IMIN,IMAX = BOUNDS ON DO-LOOP INDEX -- FIRST AND LAST C ELEMENTS OF IEND TO BE INCREMENTED C I = DO-LOOP INDEX C ITEMP = TEMPORARY STORAGE LOCATION C K = KK C C INITIALIZATION C N(1) = I1 N(2) = I2 N(3) = I3 C C SET UP NFT C DO 2 I = 1,3 N1 = N(I) INDX = MOD(I,3) + 1 N2 = N(INDX) INDX = IEND(N1) + 1 C C FIND THE INDEX OF N2 AS A NEIGHBOR OF N1 C 1 INDX = INDX - 1 IF (IADJ(INDX) .NE. N2) GO TO 1 2 NFT(I) = INDX + 1 C C ORDER THE VERTICES BY DECREASING MAGNITUDE - C N(IP(I+1)) PRECEDES N(IP(I)) IN IEND FOR C I = 1,2. C IP1 = 1 IP2 = 2 IP3 = 3 IF ( N(2) .LE. N(1) ) GO TO 3 IP1 = 2 IP2 = 1 3 IF ( N(3) .LE. N(IP1) ) GO TO 4 IP3 = IP1 IP1 = 3 4 IF ( N(IP3) .LE. N(IP2) ) GO TO 5 ITEMP = IP2 IP2 = IP3 IP3 = ITEMP C C ADD NODE K TO THE ADJACENCY LISTS OF EACH VERTEX AND C UPDATE IEND. FOR EACH VERTEX, A SET OF IADJ ELEMENTS C IS SHIFTED DOWNWARD AND K IS INSERTED. SHIFTING STARTS C AT THE END OF THE ARRAY. C 5 KM1 = K - 1 NL = IEND(KM1) NF = NFT(IP1) IF (NF .LE. NL) CALL SHIFTD(NF,NL,3, IADJ) IADJ(NF+2) = K IMIN = N(IP1) IMAX = KM1 DO 6 I = IMIN,IMAX 6 IEND(I) = IEND(I) + 3 C NL = NF - 1 NF = NFT(IP2) CALL SHIFTD(NF,NL,2, IADJ) IADJ(NF+1) = K IMAX = IMIN - 1 IMIN = N(IP2) DO 7 I = IMIN,IMAX 7 IEND(I) = IEND(I) + 2 C NL = NF - 1 NF = NFT(IP3) CALL SHIFTD(NF,NL,1, IADJ) IADJ(NF) = K IMAX = IMIN - 1 IMIN = N(IP3) DO 8 I = IMIN,IMAX 8 IEND(I) = IEND(I) + 1 C C ADD NODE K TO IEND AND ITS NEIGHBORS TO IADJ C INDX = IEND(KM1) IEND(K) = INDX + 3 DO 9 I = 1,3 INDX = INDX + 1 9 IADJ(INDX) = N(I) RETURN END SUBROUTINE INTRC0 (N,PLAT,PLON,X,Y,Z,W,IADJ,IEND, IST, . PW,IER) INTEGER N, IADJ(1), IEND(N), IST, IER REAL PLAT, PLON, X(N), Y(N), Z(N), W(N), PW C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF A SET OF NODES ON THE UNIT C SPHERE, ALONG WITH DATA VALUES AT THE NODES, THIS SUB- C ROUTINE COMPUTES THE VALUE AT A POINT P OF A CONTINUOUS C FUNCTION WHICH INTERPOLATES THE DATA VALUES. THE INTERP- C OLATORY FUNCTION IS LINEAR ON EACH UNDERLYING TRIANGLE C (PLANAR TRIANGLE WITH THE SAME VERTICES AS A SPHERICAL C TRIANGLE). IF P IS NOT CONTAINED IN A TRIANGLE, AN EX- C TRAPOLATED VALUE IS TAKEN TO BE THE INTERPOLATED VALUE AT C THE NEAREST POINT OF THE TRIANGULATION BOUNDARY. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIANGU- C LATION. N .GE. 3. C C PLAT,PLON - LATITUDE AND LONGITUDE OF P IN C RADIANS. C C X,Y,Z - VECTORS CONTAINING CARTESIAN C COORDINATES OF THE NODES. C C W - VECTOR CONTAINING DATA VALUES C AT THE NODES. W(I) IS ASSOCI- C ATED WITH (X(I),Y(I),Z(I)) FOR C I = 1,...,N. C C IADJ,IEND - TRIANGULATION DATA STRUCTURE C CREATED BY SUBROUTINE TRMESH. C C IST - INDEX OF THE STARTING NODE IN C THE SEARCH FOR A TRIANGLE CON- C TAINING P. 1 .LE. IST .LE. N. C THE OUTPUT VALUE OF IST FROM A C PREVIOUS CALL MAY BE A GOOD C CHOICE. C C INPUT PARAMETERS OTHER THAN IST ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - IST - INDEX OF ONE OF THE VERTICES OF C THE TRIANGLE CONTAINING P (OR C NEAREST P) UNLESS IER = -1 OR C IER = -2. C C PW - VALUE OF THE INTERPOLATORY C FUNCTION AT P IF IER .LE. 0. C C IER - ERROR INDICATOR C IER = 0 IF INTERPOLATION WAS C PERFORMED SUCCESSFULLY. C IER = 1 IF EXTRAPOLATION WAS C PERFORMED SUCCESSFULLY. C IER = -1 IF N .LT. 3 OR IST IS C OUT OF RANGE. C IER = -2 IF THE NODES ARE COL- C LINEAR. C IER = -3 IF P IS NOT IN A TRI- C ANGLE AND THE ANGLE BE- C TWEEN P AND THE NEAREST C BOUNDARY POINT IS AT C LEAST 90 DEGREES. C C MODULES REFERENCED BY INTRC0 - TRFIND C C INTRINSIC FUNCTIONS CALLED BY INTRC0 - COS, SIN C C*********************************************************** C INTEGER I1, I2, I3, N1, N2, INDX REAL P(3), P1(3), P2(3), P3(3), B1, B2, B3, SUM, . S12, PTN1, PTN2 C C LOCAL PARAMETERS - C C I1,I2,I3 = VERTEX INDICES RETURNED BY TRFIND C N1,N2 = ENDPOINTS OF A BOUNDARY ARC WHICH IS VISIBLE C FROM P WHEN P IS NOT CONTAINED IN A TRIANGLE C INDX = IADJ INDEX OF N1 AS A NEIGHBOR OF N2 C P = CARTESIAN COORDINATES OF P C P1,P2,P3 = CARTESIAN COORDINATES OF I1, I2, AND I3 C B1,B2,B3 = BARYCENTRIC COORDINATES OF THE CENTRAL PROJEC- C TION OF P ONTO THE UNDERLYING PLANAR TRIANGLE C OR (B1 AND B2) PROJECTION OF Q ONTO THE C UNDERLYING LINE SEGMENT N1-N2 WHEN P IS C EXTERIOR -- UNNORMALIZED COORDINATES ARE C COMPUTED BY TRFIND WHEN P IS IN A TRIANGLE C SUM = QUANTITY USED TO NORMALIZE THE BARYCENTRIC C COORDINATES C S12 = SCALAR PRODUCT (N1,N2) C PTN1 = SCALAR PRODUCT (P,N1) C PTN2 = SCALAR PRODUCT (P,N2) C IF (N .LT. 3 .OR. IST .LT. 1 .OR. IST .GT. N) . GO TO 5 C C TRANSFORM (PLAT,PLON) TO CARTESIAN COORDINATES C P(1) = COS(PLAT)*COS(PLON) P(2) = COS(PLAT)*SIN(PLON) P(3) = SIN(PLAT) C C FIND THE VERTEX INDICES OF A TRIANGLE CONTAINING P C CALL TRFIND(IST,P,X,Y,Z,IADJ,IEND, B1,B2,B3,I1,I2,I3) IF (I1 .EQ. 0) GO TO 6 IST = I1 IF (I3 .EQ. 0) GO TO 1 C C P IS CONTAINED IN THE TRIANGLE (I1,I2,I3). STORE THE C VERTEX COORDINATES IN LOCAL VARIABLES C P1(1) = X(I1) P1(2) = Y(I1) P1(3) = Z(I1) P2(1) = X(I2) P2(2) = Y(I2) P2(3) = Z(I2) P3(1) = X(I3) P3(2) = Y(I3) P3(3) = Z(I3) C C NORMALIZE THE BARYCENTRIC COORDINATES C SUM = B1 + B2 + B3 B1 = B1/SUM B2 = B2/SUM B3 = B3/SUM PW = B1*W(I1) + B2*W(I2) + B3*W(I3) IER = 0 RETURN C C P IS EXTERIOR TO THE TRIANGULATION, AND I1 AND I2 ARE C BOUNDARY NODES WHICH ARE VISIBLE FROM P. SET PW TO THE C INTERPOLATED VALUE AT Q WHERE Q IS THE CLOSEST BOUNDARY C POINT TO P. C C TRAVERSE THE BOUNDARY STARTING FROM THE RIGHTMOST VISIBLE C NODE I1. C 1 N1 = I1 PTN1 = P(1)*X(N1) + P(2)*Y(N1) + P(3)*Z(N1) IF (I1 .NE. I2) GO TO 3 C C ALL BOUNDARY NODES ARE VISIBLE FROM P. FIND A BOUNDARY C ARC N1->N2 SUCH THAT P LEFT (N2 X N1)->N1 C C COUNTERCLOCKWISE BOUNDARY TRAVERSAL -- C SET N2 TO THE FIRST NEIGHBOR OF N1 C 2 INDX = 1 IF (N1 .NE. 1) INDX = IEND(N1-1) + 1 N2 = IADJ(INDX) C C COMPUTE INNER PRODUCTS (N1,N2) AND (P,N2), AND COMPUTE C B2 = DET(P,N1,N2 X N1) C S12 = X(N1)*X(N2) + Y(N1)*Y(N2) + Z(N1)*Z(N2) PTN2 = P(1)*X(N2) + P(2)*Y(N2) + P(3)*Z(N2) B2 = PTN2 - S12*PTN1 IF (B2 .LE. 0.) GO TO 3 C C P RIGHT (N2 X N1)->N1 -- ITERATE C N1 = N2 I1 = N1 PTN1 = PTN2 GO TO 2 C C P LEFT (N2 X N1)->N1 WHERE N2 IS THE FIRST NEIGHBOR OF P1 C CLOCKWISE BOUNDARY TRAVERSAL -- C 3 N2 = N1 PTN2 = PTN1 C C SET N1 TO THE LAST NEIGHBOR OF N2 AND TEST FOR TERMINATION C INDX = IEND(N2) - 1 N1 = IADJ(INDX) IF (N1 .EQ. I1) GO TO 7 C C COMPUTE INNER PRODUCTS (N1,N2) AND (P,N1) C S12 = X(N1)*X(N2) + Y(N1)*Y(N2) + Z(N1)*Z(N2) PTN1 = P(1)*X(N1) + P(2)*Y(N1) + P(3)*Z(N1) C C COMPUTE B2 = DET(P,N1,N2 X N1) = DET(Q,N1,N2 X N1)*(P,Q) C B2 = PTN2 - S12*PTN1 IF (B2 .LE. 0.) GO TO 3 C C COMPUTE B1 = DET(P,N2 X N1,N2) = DET(Q,N2 X N1,N2)*(P,Q) C B1 = PTN1 - S12*PTN2 IF (B1 .GT. 0.) GO TO 4 C C Q = N2 C PW = W(N2) IER = 1 RETURN C C P STRICTLY LEFT (N2 X N1)->N2 AND P STRICTLY LEFT C N1->(N2 X N1). THUS Q LIES ON THE INTERIOR OF N1->N2. C NORMALIZE THE COORDINATES AND COMPUTE PW. C 4 SUM = B1 + B2 PW = (B1*W(N1) + B2*W(N2))/SUM IER = 1 RETURN C C N OR IST OUT OF RANGE C 5 IER = -1 RETURN C C COLLINEAR NODES C 6 IER = -2 RETURN C C THE ANGULAR DISTANCE BETWEEN P AND THE CLOSEST BOUNDARY C POINT TO P IS AT LEAST 90 DEGREES C 7 IER = -3 RETURN END SUBROUTINE INTRC1 (N,PLAT,PLON,X,Y,Z,W,IADJ,IEND, . IFLAG,GRAD, IST, PW,IER) INTEGER N, IADJ(1), IEND(N), IFLAG, IST, IER REAL PLAT, PLON, X(N), Y(N), Z(N), W(N), . GRAD(3,N), PW C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF A SET OF NODES ON THE UNIT C SPHERE, ALONG WITH DATA VALUES AT THE NODES, THIS SUB- C ROUTINE CONSTRUCTS THE VALUE AT A POINT P OF A ONCE CON- C TINUOUSLY DIFFERENTIABLE FUNCTION WHICH INTERPOLATES THE C DATA VALUES. IF THE TRIANGULATION DOES NOT COVER THE C ENTIRE SPHERE, THE SURFACE IS EXTENDED CONTINUOUSLY BEYOND C THE BOUNDARY ALLOWING EXTRAPOLATION. C C INPUT PARAMETERS - N - NUMBER OF NODES. N .GE. 3 AND C N .GE. 7 IF IFLAG = 0. C C PLAT,PLON - LATITUDE AND LONGITUDE OF P IN C RADIANS. C C X,Y,Z - VECTORS CONTAINING CARTESIAN C COORDINATES OF THE NODES. C C W - VECTOR CONTAINING DATA VALUES C AT THE NODES. W(I) IS ASSOCI- C ATED WITH (X(I),Y(I),Z(I)) FOR C I = 1,...,N. C C IADJ,IEND - DATA STRUCTURE REPRESENTING THE C TRIANGULATION. SEE SUBROUTINE C TRMESH. C C IFLAG - OPTION INDICATOR C IFLAG = 0 IF INTRC1 IS TO PRO- C VIDE ESTIMATED GRAD- C IENTS (FROM GRADL). C N .GE. 7 IN THIS C CASE. C IFLAG = 1 IF GRADIENTS ARE PRO- C VIDED IN GRAD. THIS C IS MORE EFFICIENT IF C INTRC1 IS TO BE C CALLED SEVERAL TIMES. C C GRAD - ARRAY DIMENSIONED 3 BY N WHOSE C I-TH COLUMN CONTAINS AN ESTI- C MATED GRADIENT AT NODE I IF C IFLAG = 1 (SEE SUBROUTINE C GRADL). GRAD MAY BE A DUMMY C VARIABLE (NOT USED) IF IFLAG C = 0. C C IST - INDEX OF THE STARTING NODE IN C THE SEARCH FOR A TRIANGLE CON- C TAINING P. 1 .LE. IST .LE. N. C THE OUTPUT VALUE OF IST FROM A C PREVIOUS CALL MAY BE A GOOD C CHOICE. C C INPUT PARAMETERS OTHER THAN IST ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - IST - INDEX OF ONE OF THE VERTICES OF C THE TRIANGLE CONTAINING P (OR C NEAREST P) UNLESS IER = -1 OR C IER = -2. C C OUTPUT PARAMETERS - PW - INTERPOLATED VALUE AT P IF C IER .GE. 0. C C IER - ERROR INDICATOR C IER = 0 IF INTERPOLATION WAS C PERFORMED SUCCESSFULLY. C IER = 1 IF EXTRAPOLATION WAS C PERFORMED SUCCESSFULLY. C IER = -1 IF N, IFLAG, OR IST IS C OUT OF RANGE. C IER = -2 IF THE NODES ARE COL- C LINEAR. C IER = -3 IF THE ANGULAR DISTANCE C BETWEEN P AND THE NEAR- C EST POINT OF THE TRIAN- C GULATION IS AT LEAST 90 C DEGREES. C C MODULES REFERENCED BY INTRC1 - TRFIND, WVAL, ARCINT, C ARCLEN, C (AND OPTIONALLY) GRADL, GETNP, CONSTR, C APLYR, SETUP, GIVENS, C ROTATE, APLYRT C C INTRINSIC FUNCTIONS CALLED BY INTRC1 - COS, SIN, SQRT C C*********************************************************** C INTEGER NN, I1, I2, I3, I, IERR, N1, N2, INDX REAL P(3), P1(3), P2(3), P3(3), W1, W2, W3, G1(3), . G2(3), G3(3), B1, B2, B3, SUM, DUM(3), S12, . PTN1, PTN2, Q(3), QNORM, WQ, GQ(3), A, PTGQ, . GQN C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C I1,I2,I3 = VERTEX INDICES RETURNED BY TRFIND C I = DO-LOOP INDEX C IERR = ERROR FLAG FOR CALLS TO GRADL C N1,N2 = INDICES OF THE ENDPOINTS OF A BOUNDARY ARC WHEN C P IS EXTERIOR (NOT CONTAINED IN A TRIANGLE) C INDX = IADJ INDEX OF N2 AS A NEIGHBOR OF N1 OR VICE- C VERSA C P = CARTESIAN COORDINATES OF P C P1,P2,P3 = CARTESIAN COORDINATES OF THE VERTICES I1, I2, C AND I3, OR (P1 AND P2) COORDINATES OF N1 AND C N2 IF P IS EXTERIOR C W1,W2,W3 = DATA VALUES ASSOCIATED WITH I1, I2, AND I3, OR C (W1 AND W2) VALUES ASSOCIATED WITH N1 AND C N2 IF P IS EXTERIOR C G1,G2,G3 = GRADIENTS AT I1, I2, AND I3, OR (G1 AND G2) AT C N1 AND N2 C B1,B2,B3 = BARYCENTRIC COORDINATES OF THE CENTRAL PROJEC- C TION OF P ONTO THE UNDERLYING PLANAR TRIANGLE C OR (B1 AND B2) PROJECTION OF Q ONTO THE C UNDERLYING LINE SEGMENT N1-N2 WHEN P IS C EXTERIOR -- UNNORMALIZED COORDINATES ARE C COMPUTED BY TRFIND WHEN P IS IN A TRIANGLE C SUM = QUANTITY USED TO NORMALIZE THE BARYCENTRIC C COORDINATES C DUM = DUMMY PARAMETER FOR WVAL AND ARCINT C S12 = SCALAR PRODUCT (N1,N2) -- FACTOR OF B1 AND B2 C PTN1 = SCALAR PRODUCT (P,N1) -- FACTOR OF B1 AND B2 C PTN2 = SCALAR PRODUCT (P,N2) -- FACTOR OF B1 AND B2 C Q = CLOSEST BOUNDARY POINT TO P WHEN P IS EXTERIOR C QNORM = FACTOR USED TO NORMALIZE Q C WQ,GQ = INTERPOLATED VALUE AND GRADIENT AT Q C A = ANGULAR SEPARATION BETWEEN P AND Q C PTGQ = SCALAR PRODUCT (P,GQ) -- FACTOR OF THE COMPONENT C OF GQ IN THE DIRECTION Q->P C GQN = NEGATIVE OF THE COMPONENT OF GQ IN THE DIRECTION C Q->P C NN = N IF (NN .LT. 3 .OR. (IFLAG .NE. 1 .AND. NN .LT. 7) . .OR. IFLAG .LT. 0 .OR. IFLAG .GT. 1 .OR. . IST .LT. 1 .OR. IST .GT. NN) GO TO 16 C C TRANSFORM (PLAT,PLON) TO CARTESIAN COORDINATES C P(1) = COS(PLAT)*COS(PLON) P(2) = COS(PLAT)*SIN(PLON) P(3) = SIN(PLAT) C C LOCATE P WITH RESPECT TO THE TRIANGULATION C CALL TRFIND(IST,P,X,Y,Z,IADJ,IEND, B1,B2,B3,I1,I2,I3) IF (I1 .EQ. 0) GO TO 17 IST = I1 IF (I3 .EQ. 0) GO TO 4 C C P IS CONTAINED IN THE TRIANGLE (I1,I2,I3). STORE THE DATA C VALUES AND VERTEX COORDINATES IN LOCAL VARIABLES. C W1 = W(I1) W2 = W(I2) W3 = W(I3) P1(1) = X(I1) P1(2) = Y(I1) P1(3) = Z(I1) P2(1) = X(I2) P2(2) = Y(I2) P2(3) = Z(I2) P3(1) = X(I3) P3(2) = Y(I3) P3(3) = Z(I3) IF (IFLAG .NE. 1) GO TO 2 C C GRADIENTS ARE USER-PROVIDED C DO 1 I = 1,3 G1(I) = GRAD(I,I1) G2(I) = GRAD(I,I2) 1 G3(I) = GRAD(I,I3) GO TO 3 C C COMPUTE GRADIENT ESTIMATES AT THE VERTICES C 2 CALL GRADL(NN,I1,X,Y,Z,W,IADJ,IEND, G1,IERR) IF (IERR .LT. 0) GO TO 17 CALL GRADL(NN,I2,X,Y,Z,W,IADJ,IEND, G2,IERR) IF (IERR .LT. 0) GO TO 17 CALL GRADL(NN,I3,X,Y,Z,W,IADJ,IEND, G3,IERR) IF (IERR .LT. 0) GO TO 17 C C NORMALIZE THE COORDINATES C 3 SUM = B1 + B2 + B3 B1 = B1/SUM B2 = B2/SUM B3 = B3/SUM CALL WVAL(B1,B2,B3,P1,P2,P3,W1,W2,W3,G1,G2,G3,0, PW, . DUM) IER = 0 RETURN C C P IS EXTERIOR TO THE TRIANGULATION, AND I1 AND I2 ARE C BOUNDARY NODES WHICH ARE VISIBLE FROM P. EXTRAPOLATE TO C P BY LINEAR (WITH RESPECT TO ARC-LENGTH) INTERPOLATION C OF THE VALUE AND DIRECTIONAL DERIVATIVE (GRADIENT COMP- C ONENT IN THE DIRECTION Q->P) OF THE INTERPOLATORY C SURFACE AT Q WHERE Q IS THE CLOSEST BOUNDARY POINT TO P. C C DETERMINE Q BY TRAVERSING THE BOUNDARY STARTING FROM I1 C 4 N1 = I1 PTN1 = P(1)*X(N1) + P(2)*Y(N1) + P(3)*Z(N1) IF (I1 .NE. I2) GO TO 6 C C ALL BOUNDARY NODES ARE VISIBLE FROM P. FIND A BOUNDARY C ARC N1->N2 SUCH THAT P LEFT (N2 X N1)->N1 C C COUNTERCLOCKWISE BOUNDARY TRAVERSAL -- C SET N2 TO THE FIRST NEIGHBOR OF N1 C 5 INDX = 1 IF (N1 .NE. 1) INDX = IEND(N1-1) + 1 N2 = IADJ(INDX) C C COMPUTE INNER PRODUCTS (N1,N2) AND (P,N2), AND COMPUTE C B2 = DET(P,N1,N2 X N1) C S12 = X(N1)*X(N2) + Y(N1)*Y(N2) + Z(N1)*Z(N2) PTN2 = P(1)*X(N2) + P(2)*Y(N2) + P(3)*Z(N2) B2 = PTN2 - S12*PTN1 IF (B2 .LE. 0.) GO TO 6 C C P RIGHT (N2 X N1)->N1 -- ITERATE C N1 = N2 I1 = N1 PTN1 = PTN2 GO TO 5 C C P LEFT (N2 X N1)->N1 WHERE N2 IS THE FIRST NEIGHBOR OF N1 C CLOCKWISE BOUNDARY TRAVERSAL -- C 6 N2 = N1 PTN2 = PTN1 C C SET N1 TO THE LAST NEIGHBOR OF N2 AND TEST FOR TERMINATION C INDX = IEND(N2) - 1 N1 = IADJ(INDX) IF (N1 .EQ. I1) GO TO 18 C C COMPUTE INNER PRODUCTS (N1,N2) AND (P,N1) C S12 = X(N1)*X(N2) + Y(N1)*Y(N2) + Z(N1)*Z(N2) PTN1 = P(1)*X(N1) + P(2)*Y(N1) + P(3)*Z(N1) C C COMPUTE B2 = DET(P,N1,N2 X N1) = DET(Q,N1,N2 X N1)*(P,Q) C B2 = PTN2 - S12*PTN1 IF (B2 .LE. 0.) GO TO 6 C C COMPUTE B1 = DET(P,N2 X N1,N2) = DET(Q,N2 X N1,N2)*(P,Q) C B1 = PTN1 - S12*PTN2 IF (B1 .GT. 0.) GO TO 10 C C Q = N2. STORE VALUE, COORDINATES, AND AND GRADIENT AT Q. C WQ = W(N2) Q(1) = X(N2) Q(2) = Y(N2) Q(3) = Z(N2) IF (IFLAG .NE. 1) GO TO 8 DO 7 I = 1,3 7 GQ(I) = GRAD(I,N2) GO TO 9 8 CALL GRADL(NN,N2,X,Y,Z,W,IADJ,IEND, GQ,IERR) IF (IERR .LT. 0) GO TO 17 C C EXTRAPOLATE TO P -- PW = WQ + A*(GQ,Q X (PXQ)/SIN(A)) C WHERE A IS THE ANGULAR SEPARATION BETWEEN Q AND P, C AND SIN(A) IS THE MAGNITUDE OF P X Q. C 9 A = ARCLEN(Q,P) PTGQ = P(1)*GQ(1) + P(2)*GQ(2) + P(3)*GQ(3) PW = WQ IF (A .NE. 0.) PW = PW + PTGQ*A/SIN(A) IER = 1 RETURN C C P STRICTLY LEFT (N2 X N1)->N2 AND P STRICTLY LEFT C N1->(N2 X N1). THUS Q LIES ON THE INTERIOR OF N1->N2. C STORE DATA VALUES AND COORDINATES OF N1 AND N2 IN C LOCAL VARIABLES C 10 W1 = W(N1) W2 = W(N2) P1(1) = X(N1) P1(2) = Y(N1) P1(3) = Z(N1) P2(1) = X(N2) P2(2) = Y(N2) P2(3) = Z(N2) C C COMPUTE THE CENTRAL PROJECTION OF Q ONTO P2-P1 AND C NORMALIZE TO OBTAIN Q C QNORM = 0. DO 11 I = 1,3 Q(I) = B1*P1(I) + B2*P2(I) 11 QNORM = QNORM + Q(I)*Q(I) QNORM = SQRT(QNORM) DO 12 I = 1,3 12 Q(I) = Q(I)/QNORM C C STORE OR COMPUTE GRADIENTS AT N1 AND N2 C IF (IFLAG .NE. 1) GO TO 14 DO 13 I = 1,3 G1(I) = GRAD(I,N1) 13 G2(I) = GRAD(I,N2) GO TO 15 14 CALL GRADL(NN,N1,X,Y,Z,W,IADJ,IEND, G1,IERR) IF (IERR .LT. 0) GO TO 17 CALL GRADL(NN,N2,X,Y,Z,W,IADJ,IEND, G2,IERR) IF (IERR .LT. 0) GO TO 17 C C COMPUTE AN INTERPOLATED VALUE AND NORMAL GRADIENT- C COMPONENT AT Q C 15 CALL ARCINT(Q,P1,P2,W1,W2,G1,G2, WQ,DUM,GQN) C C EXTRAPOLATE TO P -- THE NORMAL GRADIENT COMPONENT GQN IS C THE NEGATIVE OF THE COMPONENT IN THE DIRECTION Q->P. C PW = WQ - GQN*ARCLEN(Q,P) IER = 1 RETURN C C N, IFLAG, OR IST OUT OF RANGE C 16 IER = -1 RETURN C C COLLINEAR NODES C 17 IER = -2 RETURN C C THE DISTANCE BETWEEN P AND THE CLOSEST BOUNDARY POINT C IS AT LEAST 90 DEGREES C 18 IER = -3 RETURN END SUBROUTINE PERMUT (N,IP, A ) INTEGER N, IP(N) REAL A(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE APPLIES A SET OF PERMUTATIONS TO A VECTOR. C C INPUT PARAMETERS - N - LENGTH OF A AND IP. C C IP - VECTOR CONTAINING THE SEQUENCE OF C INTEGERS 1,...,N PERMUTED IN THE C SAME FASHION THAT A IS TO BE PER- C MUTED. C C A - VECTOR TO BE PERMUTED. C C N AND IP ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - A - REORDERED VECTOR REFLECTING THE C PERMUTATIONS DEFINED BY IP. C C MODULES REFERENCED BY PERMUT - NONE C C*********************************************************** C INTEGER NN, K, J, IPJ REAL TEMP C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C K = INDEX FOR IP AND FOR THE FIRST ELEMENT OF A IN A C PERMUTATION C J = INDEX FOR IP AND A, J .GE. K C IPJ = IP(J) C TEMP = TEMPORARY STORAGE FOR A(K) C NN = N IF (NN .LT. 2) RETURN K = 1 C C LOOP ON PERMUTATIONS C 1 J = K TEMP = A(K) C C APPLY PERMUTATION TO A. IP(J) IS MARKED (MADE NEGATIVE) C AS BEING INCLUDED IN THE PERMUTATION. C 2 IPJ = IP(J) IP(J) = -IPJ IF (IPJ .EQ. K) GO TO 3 A(J) = A(IPJ) J = IPJ GO TO 2 3 A(J) = TEMP C C SEARCH FOR AN UNMARKED ELEMENT OF IP C 4 K = K + 1 IF (K .GT. NN) GO TO 5 IF (IP(K) .GT. 0) GO TO 1 GO TO 4 C C ALL PERMUTATIONS HAVE BEEN APPLIED. UNMARK IP. C 5 DO 6 K = 1,NN 6 IP(K) = -IP(K) RETURN END SUBROUTINE PRJCT (PX,PY,PZ,OX,OY,OZ,EX,EY,EZ, . VX,VY,VZ, ISW, X,Y,IER) INTEGER ISW, IER REAL PX, PY, PZ, OX, OY, OZ, EX, EY, EZ, VX, VY, . VZ, X, Y C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE PROJECTS A POINT P ONTO THE PLANE CON- C TAINING POINT O WHOSE NORMAL IS THE LINE CONTAINING O AND C E WHERE P, O, AND E ARE POINTS IN EUCLIDEAN 3-SPACE. C C INPUT PARAMETERS - PX,PY,PZ - CARTESIAN COORDINATES OF THE C POINT P TO BE MAPPED ONTO C THE PLANE. C C OX,OY,OZ - COORDINATES OF O (THE ORIGIN C OF A COORDINATE SYSTEM IN C THE PLANE OF PROJECTION). C C EX,EY,EZ - COORDINATES OF THE EYE-POSI- C TION E DEFINING THE NORMAL C TO THE PLANE AND THE LINE OF C SIGHT FOR THE PROJECTION. E C MUST NOT COINCIDE WITH O OR C P, AND THE ANGLE BETWEEN THE C VECTORS O-E AND P-E MUST BE C LESS THAN 90 DEGREES. NOTE C THAT E AND P MAY LIE ON OP- C POSITE SIDES OF THE PLANE. C C VX,VY,VZ - COORDINATES OF A POINT V C WHICH DEFINES THE POSITIVE C Y-AXIS OF AN X-Y COORDINATE C SYSTEM IN THE PLANE AS THE C HALF-LINE CONTAINING O AND C THE PROJECTION OF O+V ONTO C THE PLANE. THE POSITIVE X- C AXIS HAS DIRECTION DEFINED C BY THE CROSS PRODUCT V X C E-O. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C ISW - SWITCH WHICH MUST BE SET TO C 1 ON THE FIRST CALL AND WHEN C THE VALUES OF OX,OY,OZ,EX, C EY,EZ,VX,VY, OR VZ HAVE BEEN C ALTERED SINCE A PREVIOUS C CALL. IF ISW .NE. 1, IT IS C ASSUMED THAT ONLY THE COORD- C INATES OF P HAVE CHANGED C SINCE A PREVIOUS CALL. PRE- C VIOUSLY STORED QUANTITIES C ARE USED FOR INCREASED EFF- C ICIENCY IN THIS CASE. C C OUTPUT PARAMETERS - ISW - RESET TO 0 IF PRJCT WAS CALLED C WITH ISW = 1, UNCHANGED OTHER- C WISE. C C X,Y - COORDINATES OF THE POINT IN THE C PLANE WHICH LIES ON THE LINE C CONTAINING E AND P. THE COORD- C INATE SYSTEM IS DEFINED BY VX, C VY, AND VZ. X AND Y ARE NOT C CHANGED IF IER .NE. 0. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF THE INNER PRODUCT OF C O-E WITH P-E IS NOT POS- C ITIVE IMPLYING THAT E IS C TOO CLOSE TO THE PLANE. C IER = 2 IF O, E, AND O+V ARE C COLLINEAR. SEE THE DES- C CRIPTION OF VX,VY,VZ. C C MODULES REFERENCED BY PRJCT - NONE C C INTRINSIC FUNCTION CALLED BY PRJCT - SQRT C C*********************************************************** C REAL XE, YE, ZE, XOE, YOE, ZOE, OES, S, XV, YV, ZV, . SC, XH, YH, ZH, XEP, YEP, ZEP, XW, YW, ZW COMMON/PROCOM/XE, YE, ZE, XOE, YOE, ZOE, OES, . XV, YV, ZV, XH, YH, ZH C C LOCAL PARAMETERS - C C XE,YE,ZE = LOCAL COPIES OF EX, EY, EZ C XOE,YOE,ZOE = COMPONENTS OF THE VECTOR OE FROM O TO E C OES = NORM SQUARED OF OE -- (OE,OE) C S = SCALE FACTOR FOR COMPUTING PROJECTIONS C XV,YV,ZV = COMPONENTS OF A UNIT VECTOR VN DEFINING THE C POSITIVE Y-AXIS IN THE PLANE C SC = SCALE FACTOR FOR NORMALIZING VN AND HN C XH,YH,ZH = COMPONENTS OF A UNIT VECTOR HN DEFINING THE C POSITIVE X-AXIS IN THE PLANE C XEP,YEP,ZEP = COMPONENTS OF THE VECTOR EP FROM E TO P C XW,YW,ZW = COMPONENTS OF THE VECTOR W FROM O TO THE C PROJECTION OF P ONTO THE PLANE C C THE COMMON BLOCK IS USED TO SAVE LOCAL PARAMETERS C BETWEEN CALLS. C IF (ISW .NE. 1) GO TO 1 ISW = 0 C C SET THE COORDINATES OF E TO LOCAL VARIABLES, COMPUTE C OE = E - O AND OES C XE = EX YE = EY ZE = EZ XOE = XE - OX YOE = YE - OY ZOE = ZE - OZ OES = XOE*XOE + YOE*YOE + ZOE*ZOE IF (OES .EQ. 0.) GO TO 2 C C COMPUTE S = (OE,V)/OES AND VN = V - S*OE C S = (XOE*VX + YOE*VY + ZOE*VZ)/OES XV = VX - S*OEX YV = VY - S*OEY ZV = VZ - S*OEZ C C NORMALIZE VN TO A UNIT VECTOR C SC = XV*XV + YV*YV + ZV*ZV IF (SC .EQ. 0.) GO TO 3 SC = 1./SQRT(SC) XV = SC*XV YV = SC*YV ZV = SC*ZV C C COMPUTE HN = VN X OE (NORMALIZED) C XH = YV*ZOE - YOE*ZV YH = XOE*ZV - XV*ZOE ZH = XV*YOE - XOE*YV SC = 1./SQRT(XH*XH + YH*YH + ZH*ZH) XH = SC*XH YH = SC*YH ZH = SC*ZH C C COMPUTE EP = P - E, S = (OE,EP)/OES, AND W = OE - EP/S C 1 XEP = PX - XE YEP = PY - YE ZEP = PZ - ZE S = (XOE*XEP + YOE*YEP + ZOE*ZEP)/OES IF (S .GE. 0.) GO TO 2 XW = XOE - XEP/S YW = YOE - YEP/S ZW = ZOE - ZEP/S C C MAP W INTO X = (W,HN), Y = (W,VN) C X = XW*XH + YW*YH + ZW*ZH Y = XW*XV + YW*YV + ZW*ZV RETURN C C (OE,EP) .GE. 0. C 2 IER = 1 RETURN C C O, E, AND O+V ARE COLLINEAR C 3 IER = 2 RETURN END SUBROUTINE QSORT (N,X, IND) INTEGER N, IND(N) REAL X(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE USES AN ORDER N*LOG(N) QUICK SORT TO C SORT THE REAL ARRAY X INTO INCREASING ORDER. THE ALGOR- C ITHM IS AS FOLLOWS. IND IS INITIALIZED TO THE ORDERED C SEQUENCE OF INDICES 1,...,N, AND ALL INTERCHANGES ARE C APPLIED TO IND. X IS DIVIDED INTO TWO PORTIONS BY PICKING C A CENTRAL ELEMENT T. THE FIRST AND LAST ELEMENTS ARE COM- C PARED WITH T, AND INTERCHANGES ARE APPLIED AS NECESSARY SO C THAT THE THREE VALUES ARE IN ASCENDING ORDER. INTER- C CHANGES ARE THEN APPLIED SO THAT ALL ELEMENTS GREATER THAN C T ARE IN THE UPPER PORTION OF THE ARRAY AND ALL ELEMENTS C LESS THAN T ARE IN THE LOWER PORTION. THE UPPER AND LOWER C INDICES OF ONE OF THE PORTIONS ARE SAVED IN LOCAL ARRAYS, C AND THE PROCESS IS REPEATED ITERATIVELY ON THE OTHER C PORTION. WHEN A PORTION IS COMPLETELY SORTED, THE PROCESS C BEGINS AGAIN BY RETRIEVING THE INDICES BOUNDING ANOTHER C UNSORTED PORTION. C C INPUT PARAMETERS - N - LENGTH OF THE ARRAY X. C C X - VECTOR OF LENGTH N TO BE SORTED. C C IND - VECTOR OF LENGTH .GE. N. C C N AND X ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - IND - SEQUENCE OF INDICES 1,...,N C PERMUTED IN THE SAME FASHION AS X C WOULD BE. THUS, THE ORDERING ON C X IS DEFINED BY Y(I) = X(IND(I)). C C MODULES REFERENCED BY QSORT - NONE C C INTRINSIC FUNCTIONS CALLED BY QSORT - IFIX, FLOAT C C*********************************************************** C C NOTE -- IU AND IL MUST BE DIMENSIONED .GE. LOG(N) WHERE C LOG HAS BASE 2. C C*********************************************************** C INTEGER IU(21), IL(21) INTEGER M, I, J, K, L, IJ, IT, ITT, INDX REAL R, T C C LOCAL PARAMETERS - C C IU,IL = TEMPORARY STORAGE FOR THE UPPER AND LOWER C INDICES OF PORTIONS OF THE ARRAY X C M = INDEX FOR IU AND IL C I,J = LOWER AND UPPER INDICES OF A PORTION OF X C K,L = INDICES IN THE RANGE I,...,J C IJ = RANDOMLY CHOSEN INDEX BETWEEN I AND J C IT,ITT = TEMPORARY STORAGE FOR INTERCHANGES IN IND C INDX = TEMPORARY INDEX FOR X C R = PSEUDO RANDOM NUMBER FOR GENERATING IJ C T = CENTRAL ELEMENT OF X C IF (N .LE. 0) RETURN C C INITIALIZE IND, M, I, J, AND R C DO 1 I = 1,N 1 IND(I) = I M = 1 I = 1 J = N R = .375 C C TOP OF LOOP C 2 IF (I .GE. J) GO TO 10 IF (R .GT. .5898437) GO TO 3 R = R + .0390625 GO TO 4 3 R = R - .21875 C C INITIALIZE K C 4 K = I C C SELECT A CENTRAL ELEMENT OF X AND SAVE IT IN T C IJ = I + IFIX(R*FLOAT(J-I)) IT = IND(IJ) T = X(IT) C C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T, C INTERCHANGE IT WITH T C INDX = IND(I) IF (X(INDX) .LE. T) GO TO 5 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) C C INITIALIZE L C 5 L = J C C IF THE LAST ELEMENT OF THE ARRAY IS LESS THAN T, C INTERCHANGE IT WITH T C INDX = IND(J) IF (X(INDX) .GE. T) GO TO 7 IND(IJ) = INDX IND(J) = IT IT = INDX T = X(IT) C C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T, C INTERCHANGE IT WITH T C INDX = IND(I) IF (X(INDX) .LE. T) GO TO 7 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) GO TO 7 C C INTERCHANGE ELEMENTS K AND L C 6 ITT = IND(L) IND(L) = IND(K) IND(K) = ITT C C FIND AN ELEMENT IN THE UPPER PART OF THE ARRAY WHICH IS C NOT LARGER THAN T C 7 L = L - 1 INDX = IND(L) IF (X(INDX) .GT. T) GO TO 7 C C FIND AN ELEMENT IN THE LOWER PART OF THE ARRAY WHCIH IS C NOT SMALLER THAN T C 8 K = K + 1 INDX = IND(K) IF (X(INDX) .LT. T) GO TO 8 C C IF K .LE. L, INTERCHANGE ELEMENTS K AND L C IF (K .LE. L) GO TO 6 C C SAVE THE UPPER AND LOWER SUBSCRIPTS OF THE PORTION OF THE C ARRAY YET TO BE SORTED C IF (L-I .LE. J-K) GO TO 9 IL(M) = I IU(M) = L I = K M = M + 1 GO TO 11 C 9 IL(M) = K IU(M) = J J = L M = M + 1 GO TO 11 C C BEGIN AGAIN ON ANOTHER UNSORTED PORTION OF THE ARRAY C 10 M = M - 1 IF (M .EQ. 0) RETURN I = IL(M) J = IU(M) C 11 IF (J-I .GE. 11) GO TO 4 IF (I .EQ. 1) GO TO 2 I = I - 1 C C SORT ELEMENTS I+1,...,J. NOTE THAT 1 .LE. I .LT. J AND C J-I .LT. 11. C 12 I = I + 1 IF (I .EQ. J) GO TO 10 INDX = IND(I+1) T = X(INDX) IT = INDX INDX = IND(I) IF (X(INDX) .LE. T) GO TO 12 K = I C 13 IND(K+1) = IND(K) K = K - 1 INDX = IND(K) IF (T .LT. X(INDX)) GO TO 13 IND(K+1) = IT GO TO 12 END SUBROUTINE REORDR (N,IFLAG, A,B,C,D, IND) INTEGER N, IFLAG, IND(N) REAL A(N), B(N), C(N), D(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE USES A QUICK SORT TO REORDER THE REAL C ARRAY A INTO INCREASING ORDER. A RECORD OF THE PERMUTA- C TIONS APPLIED TO A IS STORED IN IND, AND THESE PERMUTA- C TIONS MAY BE APPLIED TO ONE, TWO, OR THREE ADDITIONAL C VECTORS BY THIS ROUTINE. ANY OTHER VECTOR V MAY BE PER- C MUTED IN THE SAME FASHION BY CALLING SUBROUTINE PERMUT C WITH N, IND, AND V AS PARAMETERS. C A SET OF NODES (X(I),Y(I),Z(I)) AND DATA VALUES W(I) C MAY BE PREPROCESSED BY REORDR FOR INCREASED EFFICIENCY IN C THE TRIANGULATION ROUTINE TRMESH. THIS PREPROCESSING IS C ALSO USEFUL FOR DETECTING DUPLICATE NODES. NOTE THAT THE C FIRST THREE NODES MUST NOT BE COLLINEAR ON INPUT TO C TRMESH. EITHER X, Y, OR Z MAY BE USED AS THE SORT KEY C (ASSOCIATED WITH A). IF THE NODAL COORDINATES ARE IN C TERMS OF LATITUDE AND LONGITUDE, IT IS USUALLY ADVAN- C TAGEOUS TO SORT ON LONGITUDE WITH THE SAME INTERCHANGES C APPLIED TO LATITUDE AND THE DATA VALUES. D WOULD BE A C DUMMY PARAMETER IN THIS CASE. C C INPUT PARAMETERS - N - NUMBER OF NODES. C C IFLAG - NUMBER OF VECTORS TO BE PER- C MUTED. C IFLAG .LE. 0 IF A, B, C, AND C D ARE TO REMAIN UNALTERED. C IFLAG = 1 IF ONLY A IS TO BE C PERMUTED. C IFLAG = 2 IF A AND B ARE TO BE C PERMUTED. C IFLAG = 3 IF A, B, AND C ARE C TO BE PERMUTED. C IFLAG .GE. 4 IF A, B, C, AND D C ARE TO BE PERMUTED. C C A,B,C,D - VECTORS OF LENGTH N TO BE C SORTED (ON THE COMPONENTS OF A) C OR DUMMY PARAMETERS, DEPENDING C ON IFLAG. C C IND - VECTOR OF LENGTH .GE. N. C C N, IFLAG, AND ANY DUMMY PARAMETERS ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - A,B,C,D - SORTED OR UNALTERED VECTORS, C DEPENDING ON IFLAG. C C IND - SEQUENCE OF INDICES 1,...,N C PERMUTED IN THE SAME FASHION C AS THE REAL VECTORS. THUS C THE ORDERING MAY BE APPLIED C TO A VECTOR V AND STORED IN C W BY SETTING W(I) = C V(IND(I)), OR V MAY BE OVER- C WRITTEN WITH THE ORDERING BY C A CALL TO PERMUT. C C MODULES REFERENCED BY REORDR - QSORT, PERMUT C C*********************************************************** C INTEGER NN, NV C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C NV = LOCAL COPY OF IFLAG C NN = N NV = IFLAG CALL QSORT (NN,A, IND) IF (NV .LE. 0) RETURN CALL PERMUT (NN,IND, A ) IF (NV .EQ. 1) RETURN CALL PERMUT (NN,IND, B ) IF (NV .EQ. 2) RETURN CALL PERMUT (NN,IND, C ) IF (NV .EQ. 3) RETURN CALL PERMUT (NN,IND, D ) RETURN END SUBROUTINE ROTATE (N,C,S, X,Y ) INTEGER N REAL C, S, X(N), Y(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C ( C S) C THIS ROUTINE APPLIES THE GIVENS ROTATION ( ) TO THE C (-S C) C (X(1) ... X(N)) C 2 BY N MATRIX ( ). THIS ROUTINE WAS TAKEN C (Y(1) ... Y(N)) C FROM LINPACK. C C INPUT PARAMETERS - N - NUMBER OF COLUMNS TO BE ROTATED. C C C,S - ELEMENTS OF THE GIVENS ROTATION. C THESE MAY BE DETERMINED BY C SUBROUTINE GIVENS. C C X,Y - VECTORS OF LENGTH .GE. N C CONTAINING THE 2-VECTORS TO BE C ROTATED. C C THE PARAMETERS N, C, AND S ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - X,Y - ROTATED VECTORS C C MODULES REFERENCED BY ROTATE - NONE C C*********************************************************** C INTEGER I REAL XI, YI C C LOCAL PARAMETERS - C C I = DO-LOOP INDEX C XI,YI = X(I), Y(I) C IF (N .LE. 0 .OR. (C .EQ. 1. .AND. S .EQ. 0.)) RETURN DO 1 I = 1,N XI = X(I) YI = Y(I) X(I) = C*XI + S*YI 1 Y(I) = -S*XI + C*YI RETURN END SUBROUTINE SETUP (XI,YI,WI,WK,S1,S2,WT, ROW) REAL XI, YI, WI, WK, S1, S2, WT, ROW(6) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE SETS UP THE I-TH ROW OF AN AUGMENTED C REGRESSION MATRIX FOR A WEIGHTED LEAST SQUARES FIT OF A C QUADRATIC FUNCTION Q(X,Y) TO A SET OF DATA VALUES WI C WHERE Q(0,0) = WK. THE FIRST 3 COLUMNS (QUADRATIC TERMS) C ARE SCALED BY 1/S2 AND THE FOURTH AND FIFTH COLUMNS (LIN- C EAR TERMS) ARE SCALED BY 1/S1. C C INPUT PARAMETERS - XI,YI - COORDINATES OF NODE I. C C WI - DATA VALUE AT NODE I. C C WK - DATA VALUE INTERPOLATED BY Q AT C THE ORIGIN. C C S1,S2 - INVERSE SCALE FACTORS. C C WT - WEIGHT FACTOR CORRESPONDING TO C THE I-TH EQUATION. C C ROW - VECTOR OF LENGTH 6. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - ROW - VECTOR CONTAINING A ROW OF THE C AUGMENTED REGRESSION MATRIX. C C MODULES REFERENCED BY SETUP - NONE C C*********************************************************** C REAL W1, W2 C C LOCAL PARAMETERS - C C W1 = WEIGHTED SCALE FACTOR FOR THE LINEAR TERMS C W2 = WEIGHTED SCALE FACTOR FOR THE QUADRATIC TERMS C W1 = WT/S1 W2 = WT/S2 ROW(1) = XI*XI*W2 ROW(2) = XI*YI*W2 ROW(3) = YI*YI*W2 ROW(4) = XI*W1 ROW(5) = YI*W1 ROW(6) = (WI-WK)*WT RETURN END SUBROUTINE SHIFTD (NFRST,NLAST,KK, IARR ) INTEGER NFRST, NLAST, KK, IARR(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE SHIFTS A SET OF CONTIGUOUS ELEMENTS OF AN C INTEGER ARRAY KK POSITIONS DOWNWARD (UPWARD IF KK .LT. 0). C THE LOOPS ARE UNROLLED IN ORDER TO INCREASE EFFICIENCY. C C INPUT PARAMETERS - NFRST,NLAST - BOUNDS ON THE PORTION OF C IARR TO BE SHIFTED. ALL C ELEMENTS BETWEEN AND C INCLUDING THE BOUNDS ARE C SHIFTED UNLESS NFRST .GT. C NLAST, IN WHICH CASE NO C SHIFT OCCURS. C C KK - NUMBER OF POSITIONS EACH C ELEMENT IS TO BE SHIFTED. C IF KK .LT. 0 SHIFT UP. C IF KK .GT. 0 SHIFT DOWN. C C IARR - INTEGER ARRAY OF LENGTH C .GE. NLAST + MAX(KK,0). C C NFRST, NLAST, AND KK ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - IARR - SHIFTED ARRAY. C C MODULES REFERENCED BY SHIFTD - NONE C C*********************************************************** C INTEGER INC, K, NF, NL, NLP1, NS, NSL, I, IBAK, INDX, . IMAX DATA INC/5/ C C LOCAL PARAMETERS - C C INC = DO-LOOP INCREMENT (UNROLLING FACTOR) -- IF INC IS C CHANGED, STATEMENTS MUST BE ADDED TO OR DELETED C FROM THE DO-LOOPS C K = LOCAL COPY OF KK C NF = LOCAL COPY OF NFRST C NL = LOCAL COPY OF NLAST C NLP1 = NL + 1 C NS = NUMBER OF SHIFTS C NSL = NUMBER OF SHIFTS DONE IN UNROLLED DO-LOOP (MULTIPLE C OF INC) C I = DO-LOOP INDEX AND INDEX FOR IARR C IBAK = INDEX FOR DOWNWARD SHIFT OF IARR C INDX = INDEX FOR IARR C IMAX = BOUND ON DO-LOOP INDEX C K = KK NF = NFRST NL = NLAST IF (NF .GT. NL .OR. K .EQ. 0) RETURN NLP1 = NL + 1 NS = NLP1 - NF NSL = INC*(NS/INC) IF ( K .LT. 0) GO TO 4 C C SHIFT DOWNWARD STARTING FROM THE BOTTOM C IF (NSL .LE. 0) GO TO 2 DO 1 I = 1,NSL,INC IBAK = NLP1 - I INDX = IBAK + K IARR(INDX) = IARR(IBAK) IARR(INDX-1) = IARR(IBAK-1) IARR(INDX-2) = IARR(IBAK-2) IARR(INDX-3) = IARR(IBAK-3) IARR(INDX-4) = IARR(IBAK-4) 1 CONTINUE C C PERFORM THE REMAINING NS-NSL SHIFTS ONE AT A TIME C 2 IBAK = NLP1 - NSL 3 IF (IBAK .LE. NF) RETURN IBAK = IBAK - 1 INDX = IBAK + K IARR(INDX) = IARR(IBAK) GO TO 3 C C SHIFT UPWARD STARTING FROM THE TOP C 4 IF (NSL .LE. 0) GO TO 6 IMAX = NLP1 - INC DO 5 I = NF,IMAX,INC INDX = I + K IARR(INDX) = IARR(I) IARR(INDX+1) = IARR(I+1) IARR(INDX+2) = IARR(I+2) IARR(INDX+3) = IARR(I+3) IARR(INDX+4) = IARR(I+4) 5 CONTINUE C C PERFORM THE REMAINING NS-NSL SHIFTS ONE AT A TIME C 6 I = NSL + NF 7 IF (I .GT. NL) RETURN INDX = I + K IARR(INDX) = IARR(I) I = I + 1 GO TO 7 END SUBROUTINE SWAP (NIN1,NIN2,NOUT1,NOUT2, IADJ,IEND ) INTEGER NIN1, NIN2, NOUT1, NOUT2, IADJ(1), IEND(1) EXTERNAL INDEX C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE SWAPS THE DIAGONALS IN A CONVEX QUADRI- C LATERAL. C C INPUT PARAMETERS - NIN1,NIN2,NOUT1,NOUT2 - NODAL INDICES C OF A PAIR OF ADJACENT TRIANGLES C WHICH FORM A CONVEX QUADRILAT- C ERAL. NOUT1 AND NOUT2 ARE CON- C NECTED BY AN ARC WHICH IS TO BE C REPLACED BY THE ARC NIN1-NIN2. C (NIN1,NOUT1,NOUT2) MUST BE TRI- C ANGLE VERTICES IN COUNTERCLOCK- C WISE ORDER. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C IADJ,IEND - TRIANGULATION DATA STRUCTURE C (SEE SUBROUTINE TRMESH). C C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ARC C REPLACEMENT. C C MODULES REFERENCED BY SWAP - INDEX, SHIFTD C C*********************************************************** C INTEGER IN(2), IO(2), IP1, IP2, J, K, NF, NL, I, . IMIN, IMAX C C LOCAL PARAMETERS - C C IN = NIN1 AND NIN2 ORDERED BY INCREASING MAGNITUDE C (THE NEIGHBORS OF IN(1) PRECEDE THOSE OF C IN(2) IN IADJ) C IO = NOUT1 AND NOUT2 IN INCREASING ORDER C IP1,IP2 = PERMUTATION OF (1,2) SUCH THAT IO(IP1) C PRECEDES IO(IP2) AS A NEIGHBOR OF IN(1) C J,K = PERMUTATION OF (1,2) USED AS INDICES OF IN C AND IO C NF,NL = IADJ INDICES BOUNDARY A PORTION OF THE ARRAY C TO BE SHIFTED C I = IEND INDEX C IMIN,IMAX = BOUNDS ON THE PORTION OF IEND TO BE INCRE- C MENTED OR DECREMENTED C IN(1) = NIN1 IN(2) = NIN2 IO(1) = NOUT1 IO(2) = NOUT2 IP1 = 1 C C ORDER THE INDICES SO THAT IN(1) .LT. IN(2) AND IO(1) .LT. C IO(2), AND CHOOSE IP1 AND IP2 SUCH THAT (IN(1),IO(IP1), C IO(IP2)) FORMS A TRIANGLE. C IF (IN(1) .LT. IN(2)) GO TO 1 IN(1) = IN(2) IN(2) = NIN1 IP1 = 2 1 IF (IO(1) .LT. IO(2)) GO TO 2 IO(1) = IO(2) IO(2) = NOUT1 IP1 = 3 - IP1 2 IP2 = 3 - IP1 IF (IO(2) .LT. IN(1)) GO TO 8 IF (IN(2) .LT. IO(1)) GO TO 12 C C IN(1) AND IO(1) PRECEDE IN(2) AND IO(2). FOR (J,K) = C (1,2) AND (2,1), DELETE IO(K) AS A NEIGHBOR OF IO(J) C BY SHIFTING A PORTION OF IADJ EITHER UP OR DOWN AND C AND INSERT IN(K) AS A NEIGHBOR OF IN(J). C DO 7 J = 1,2 K = 3 - J IF (IN(J) .GT. IO(J)) GO TO 4 C C THE NEIGHBORS OF IN(J) PRECEDE THOSE OF IO(J) -- SHIFT C DOWN BY 1 C NF = 1 + INDEX(IN(J),IO(IP1),IADJ,IEND) NL = -1 + INDEX(IO(J),IO(K),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ ) IADJ(NF) = IN(K) IMIN = IN(J) IMAX = IO(J)-1 DO 3 I = IMIN,IMAX 3 IEND(I) = IEND(I) + 1 GO TO 6 C C THE NEIGHBORS OF IO(J) PRECEDE THOSE OF IN(J) -- SHIFT C UP BY 1 C 4 NF = 1 + INDEX(IO(J),IO(K),IADJ,IEND) NL = -1 + INDEX(IN(J),IO(IP2),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ ) IADJ(NL) = IN(K) IMIN = IO(J) IMAX = IN(J) - 1 DO 5 I = IMIN,IMAX 5 IEND(I) = IEND(I) - 1 C C REVERSE (IP1,IP2) FOR (J,K) = (2,1) C 6 IP1 = IP2 IP2 = 3 - IP1 7 CONTINUE RETURN C C THE VERTICES ARE ORDERED (IO(1),IO(2),IN(1),IN(2)). C DELETE IO(2) BY SHIFTING UP BY 1 C 8 NF = 1 + INDEX(IO(1),IO(2),IADJ,IEND) NL = -1 + INDEX(IO(2),IO(1),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ ) IMIN = IO(1) IMAX = IO(2)-1 DO 9 I = IMIN,IMAX 9 IEND(I) = IEND(I) - 1 C C DELETE IO(1) BY SHIFTING UP BY 2 AND INSERT IN(2) C NF = NL + 2 NL = -1 + INDEX(IN(1),IO(IP2),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,-2, IADJ ) IADJ(NL-1) = IN(2) IMIN = IO(2) IMAX = IN(1)-1 DO 10 I = IMIN,IMAX 10 IEND(I) = IEND(I) - 2 C C SHIFT UP BY 1 AND INSERT IN(1) C NF = NL + 1 NL = -1 + INDEX(IN(2),IO(IP1),IADJ,IEND) CALL SHIFTD(NF,NL,-1, IADJ ) IADJ(NL) = IN(1) IMIN = IN(1) IMAX = IN(2)-1 DO 11 I = IMIN,IMAX 11 IEND(I) = IEND(I) - 1 RETURN C C THE VERTICES ARE ORDERED (IN(1),IN(2),IO(1),IO(2)). C DELETE IO(1) BY SHIFTING DOWN BY 1 C 12 NF = 1 + INDEX(IO(1),IO(2),IADJ,IEND) NL = -1 + INDEX(IO(2),IO(1),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ ) IMIN = IO(1) IMAX = IO(2) - 1 DO 13 I = IMIN,IMAX 13 IEND(I) = IEND(I) + 1 C C DELETE IO(2) BY SHIFTING DOWN BY 2 AND INSERT IN(1) C NL = NF - 2 NF = 1 + INDEX(IN(2),IO(IP2),IADJ,IEND) IF (NF .LE. NL) CALL SHIFTD(NF,NL,2, IADJ ) IADJ(NF+1) = IN(1) IMIN = IN(2) IMAX = IO(1) - 1 DO 14 I = IMIN,IMAX 14 IEND(I) = IEND(I) + 2 C C SHIFT DOWN BY 1 AND INSERT IN(2) C NL = NF - 1 NF = 1 + INDEX(IN(1),IO(IP1),IADJ,IEND) CALL SHIFTD(NF,NL,1, IADJ ) IADJ(NF) = IN(2) IMIN = IN(1) IMAX = IN(2) - 1 DO 15 I = IMIN,IMAX 15 IEND(I) = IEND(I) + 1 RETURN END LOGICAL FUNCTION SWPTST (N1,N2,N3,N4,X,Y,Z) INTEGER N1, N2, N3, N4 REAL X(1), Y(1), Z(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 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 POSITIVE (SWPTST = .TRUE.) IF AND C ONLY IF N4 LIES ABOVE THE PLANE (IN THE HALF-SPACE NOT C CONTAINING THE ORIGIN) DEFINED BY (N1,N2,N3), OR EQUIV- C ALENTLY, IF THE PROJECTION OF N4 ONTO THIS PLANE IS C INTERIOR TO THE CIRCUMCIRCLE OF (N1,N2,N3). THE DECISION C WILL BE NEGATIVE IF THE QUADRILATERAL IS NOT STRICTLY C CONVEX. C C INPUT PARAMETERS - N1,N2,N3,N4 - INDICES OF THE FOUR NODES C DEFINING THE QUADRILATERAL. N1 AND C N2 ARE CURRENTLY CONNECTED BY A DIAG- C ONAL ARC WHICH SHOULD BE REPLACED BY C AN ARC CONNECTING N3 TO N4 IF THE C DECISION IS MADE TO SWAP. N1, N2, N3 C MUST BE IN COUNTERCLOCKWISE ORDER. C C X,Y,Z - VECTORS OF NODAL COORDINATES. (X(I), C Y(I),Z(I)) ARE THE COORDINATES OF C NODE I FOR I = N1, N2, N3, AND N4. C C NONE OF THE INPUT PARAMETERS ARE ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - SWPTST - .TRUE. IFF THE ARC CONNECTING C N1 AND N2 IS TO BE REPLACED. C C MODULES REFERENCED BY SWPTST - NONE C C*********************************************************** C REAL X4, Y4, Z4, DX1, DX2, DX3, DY1, DY2, DY3, . DZ1, DZ2, DZ3 C C LOCAL PARAMETERS - C C X4,Y4,Z4 = COORDINATES OF N4 C DX1,DY1,DZ1 = COORDINATES OF N1 - N4 C DX2,DY2,DZ2 = COORDINATES OF N2 - N4 C DX3,DY3,DZ3 = COORDINATES OF N3 - 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) .GT. 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 ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 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 INPUT PARAMETERS - N - NUMBER OF NODES (POINTS ON THE C UNIT SPHERE) WHOSE COORDINATES C ARE TO BE TRANSFORMED. C C RLAT - N-VECTOR CONTAINING LATITUDINAL C COORDINATES OF THE NODES IN C RADIANS. C C RLON - N-VECTOR CONTAINING LONGITUDINAL C COORDINATES OF THE NODES IN C RADIANS. C C X,Y,Z - VECTORS OF LENGTH .GE. N. C C N, RLAT, AND RLON ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - X,Y,Z - CARTESIAN COORDINATES IN THE C RANGE (-1,1). X(I)**2 + C Y(I)**2 + Z(I)**2 = 1 FOR C I = 1,...,N. C C MODULES REFERENCED BY TRANS - NONE C C INTRINSIC FUNCTIONS CALLED BY TRANS - COS, SIN C C*********************************************************** C INTEGER NN, I REAL PHI, THETA, COSPHI C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C I = DO-LOOP INDEX C PHI = LATITUDE C THETA = LONGITUDE C COSPHI = COS(PHI) 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) 1 Z(I) = SIN(PHI) RETURN END SUBROUTINE TRFIND (NST,P,X,Y,Z,IADJ,IEND, B1,B2,B3, . I1,I2,I3) INTEGER NST, IADJ(1), IEND(1), I1, I2, I3 REAL P(3), X(1), Y(1), Z(1), B1, B2, B3 C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE LOCATES A POINT P RELATIVE TO A TRIANGU- C LATION OF THE CONVEX HULL OF A SET OF NODES (POINTS ON THE C UNIT SPHERE). C C INPUT PARAMETERS - NST - INDEX OF NODE AT WHICH TRFIND C BEGINS SEARCH. SEARCH TIME C DEPENDS ON THE PROXIMITY OF C NST TO P. C C P - X-, Y-, AND Z-COORDINATES (IN C THAT ORDER) OF THE POINT TO BE C LOCATED. C C X,Y,Z - VECTORS CONTAINING CARTESIAN C COORDINATES OF THE NODES IN C THE MESH. (X(I),Y(I),Z(I)) C DEFINES NODE I FOR I = 1,...,N C WHERE N .GE. 3. C C IADJ - SET OF ADJACENCY LISTS OF C NODES IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IADJ AND IEND MAY BE CREATED BY TRMESH. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - B1,B2,B3 - UNNORMALIZED BARYCENTRIC C COORDINATES OF THE CENTRAL C PROJECTION OF P ONTO THE C UNDERLYING PLANAR TRIANGLE C IF P IS IN THE CONVEX HULL C OF THE NODES, UNCHANGED IF C I1 = 0. C C I1,I2,I3 - COUNTERCLOCKWISE-ORDERED C VERTEX INDICES OF A TRI- C ANGLE CONTAINING P IF P IS C AN INTERIOR POINT. IF P IS C NOT CONTAINED IN THE CONVEX C HULL OF THE NODES, I1 AND C I2 ARE THE RIGHTMOST AND C LEFTMOST NODES WHICH ARE C VISIBLE FROM P, AND I3 = 0. C IF ALL BOUNDARY NODES ARE C VISIBLE FROM P, I1 AND I2 C COINCIDE. IF P AND ALL OF C THE NODES LIE ON A SINGLE C GREAT CIRCLE THEN I1 = I2 C = I3 = 0. C C MODULES REFERENCED BY TRFIND - NONE C C INTRINSIC FUNCTIONS CALLED BY TRFIND - MAX0, ABS C C*********************************************************** C INTEGER N0, N1, N2, N3, N4, INDX, IND, NF, NL, N1S, . N2S, NEXT, NRST, NRMAX, LUN REAL XP, YP, ZP, B3P1, S12, PTN1, PTN2, Q(3) DATA NRMAX/5/, LUN/6/ C C LOCAL PARAMETERS - C 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 N3,N4 = NODES OPPOSITE N1->N2 AND N2->N1, RESPECTIVELY C INDX,IND = INDICES FOR IADJ C NF,NL = FIRST AND LAST NEIGHBORS OF N0 IN IADJ, OR C FIRST (RIGHTMOST) AND LAST (LEFTMOST) NODES C VISIBLE FROM P WHEN P IS EXTERIOR TO THE C BOUNDARY C N1S,N2S = INITIALLY-DETERMINED VALUES OF N1 AND N2 WHEN C P IS EXTERIOR C NEXT = CANDIDATE FOR I1 OR I2 WHEN P IS EXTERIOR C NRST = NUMBER OF RESTARTS WITH NEW N0 C NRMAX = MAXIMUM ALLOWABLE NUMBER OF RESTARTS BEFORE C TERMINATING WITH AN ERROR MESSAGE C LUN = LOGICAL UNIT FOR ERROR MESSAGES C XP,YP,ZP = LOCAL VARIABLES CONTAINING P(1), P(2), AND P(3) C B3P1 = B3 + 1 -- B3P1 = 1 IFF B3 = 0 TO WITHIN THE C MACHINE PRECISION C S12 = SCALAR PRODUCT 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 DET = STATEMENT FUNCTION WHICH COMPUTES A DETERMI- C NANT -- DET(X1,...,Z0) .GE. 0 IFF (X0,Y0,Z0) C IS IN THE (CLOSED) LEFT HEMISPHERE DEFINED BY C THE PLANE CONTAINING (0,0,0), (X1,Y1,Z1), C AND (X2,Y2,Z2) WHERE LEFT IS DEFINED RELA- C TIVE TO AN OBSERVER AT (X1,Y1,Z1) FACING C (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) C C INITIALIZE VARIABLES C XP = P(1) YP = P(2) ZP = P(3) NRST = 0 N0 = MAX0(NST,1) C C FIND A CONE WITH VERTEX N0 CONTAINING P C 1 INDX = IEND(N0) NL = IADJ(INDX) INDX = 1 IF (N0 .NE. 1) INDX = IEND(N0-1) + 1 NF = IADJ(INDX) N1 = NF IF (NL .NE. 0) GO TO 3 C C N0 IS A BOUNDARY NODE. SET NL TO THE LAST NONZERO C NEIGHBOR OF N0. C IND = IEND(N0) - 1 NL = IADJ(IND) IF ( DET(X(N0),Y(N0),Z(N0),X(NF),Y(NF),Z(NF), . XP,YP,ZP) .GE. 0. ) GO TO 2 C C P IS TO THE RIGHT OF THE BOUNDARY EDGE N0->NF C N1 = N0 N2 = NF GO TO 16 2 IF ( DET(X(NL),Y(NL),Z(NL),X(N0),Y(N0),Z(N0), . XP,YP,ZP) .GE. 0. ) GO TO 4 C C P IS TO THE RIGHT OF THE BOUNDARY EDGE NL->N0 C N1 = NL N2 = N0 GO TO 16 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) .GE. 0. ) GO TO 4 INDX = INDX + 1 N1 = IADJ(INDX) IF (N1 .EQ. NL) GO TO 7 GO TO 3 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 INDX = INDX + 1 N2 = IADJ(INDX) IF ( DET(X(N0),Y(N0),Z(N0),X(N2),Y(N2),Z(N2), . XP,YP,ZP) .LT. 0. ) GO TO 8 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 7 C C P IS LEFT OF OR ON ARCS N0->NB FOR ALL NEIGHBORS NB C OF N0. TEST FOR P = +/-N0. C IF (ABS(X(N0)*XP + Y(N0)*YP + Z(N0)*ZP) .GE. 1.) . GO TO 6 C C ALL POINTS ARE COLLINEAR IFF P IS LEFT OF NB->N0 FOR C ALL NEIGHBORS NB OF N0. SEARCH THE NEIGHBORS OF N0 C IN REVERSE ORDER. NOTE -- N1 = NL AND INDX POINTS TO C NL. C 5 IF ( DET(X(N1),Y(N1),Z(N1),X(N0),Y(N0),Z(N0), . XP,YP,ZP) .LT. 0. ) GO TO 6 IF (N1 .EQ. NF) GO TO 24 INDX = INDX - 1 N1 = IADJ(INDX) GO TO 5 C C P IS TO THE RIGHT OF N1->N0, OR P = +/-N0. SET N0 TO N1 C AND START OVER. C 6 N0 = N1 NRST = NRST + 1 IF (NRST .EQ. NRMAX) GO TO 25 GO TO 1 C C P IS BETWEEN ARCS N0->N1 AND N0->NF C 7 N2 = NF C C P IS CONTAINED IN A CONE DEFINED BY LINE SEGMENTS N0-N1 C AND N0-N2 WHERE N1 IS ADJACENT TO N2 C 8 N3 = N0 9 B3 = DET(X(N1),Y(N1),Z(N1),X(N2),Y(N2),Z(N2),XP,YP,ZP) IF (B3 .GE. 0.) GO TO 13 C C SET N4 TO THE FIRST NEIGHBOR OF N2 FOLLOWING N1 C INDX = IEND(N2) IF (IADJ(INDX) .NE. N1) GO TO 10 C C N1 IS THE LAST NEIGHBOR OF N2. C SET N4 TO THE FIRST NEIGHBOR. C INDX = 1 IF (N2 .NE. 1) INDX = IEND(N2-1) + 1 N4 = IADJ(INDX) GO TO 11 C C N1 IS NOT THE LAST NEIGHBOR OF N2 C 10 INDX = INDX-1 IF (IADJ(INDX) .NE. N1) GO TO 10 N4 = IADJ(INDX+1) IF (N4 .EQ. 0) GO TO 16 C C DEFINE A NEW ARC N1->N2 WHICH INTERSECTS THE LINE C SEGMENT N0-P C 11 IF ( DET(X(N0),Y(N0),Z(N0),X(N4),Y(N4),Z(N4), . XP,YP,ZP) .GE. 0. ) GO TO 12 N3 = N2 N2 = N4 GO TO 9 12 N3 = N1 N1 = N4 GO TO 9 C C P IS IN (N1,N2,N3) UNLESS N0, N1, N2, AND P ARE COLLINEAR C OR P IS CLOSE TO -N0. C 13 B3P1 = B3 + 1. IF (B3P1 .LE. 1.) GO TO 14 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 .GE. 0. .AND. B2 .GE. 0.) GO TO 15 C C RESTART WITH N0 = N3 C NRST = NRST + 1 IF (NRST .EQ. NRMAX) GO TO 25 N0 = N3 GO TO 1 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 14 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 .GE. 0. .AND. B2 .GE. 0.) GO TO 15 C C RESTART WITH N0 = 1 C NRST = NRST + 1 IF (NRST .EQ. NRMAX) GO TO 25 N0 = 1 GO TO 1 C C P IS IN (N1,N2,N3) C 15 I1 = N1 I2 = N2 I3 = N3 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 16 N1S = N1 N2S = N2 NL = 0 C C COUNTERCLOCKWISE BOUNDARY TRAVERSAL C 17 INDX = 1 IF (N2 .NE. 1) INDX = IEND(N2-1) + 1 NEXT = IADJ(INDX) IF ( DET(X(N2),Y(N2),Z(N2),X(NEXT),Y(NEXT),Z(NEXT), . XP,YP,ZP) .LT. 0. ) GO TO 18 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 19 IF (X(NEXT)*Q(1) + Y(NEXT)*Q(2) + Z(NEXT)*Q(3) . .GE. 0.) GO TO 19 C C N1, N2, NEXT, AND P ARE NEARLY COLLINEAR, AND N2 IS C THE LEFTMOST VISIBLE NODE C NL = N2 18 N1 = N2 N2 = NEXT IF (N2 .NE. N1S) GO TO 17 C C ALL BOUNDARY NODES ARE VISIBLE C I1 = N1S I2 = N1S I3 = 0 RETURN C C N2 IS THE RIGHTMOST VISIBLE NODE C 19 NF = N2 IF (NL .NE. 0) GO TO 23 C C RESTORE INITIAL VALUES OF N1 AND N2, AND BEGIN SEARCH C FOR THE LEFTMOST VISIBLE NODE C N2 = N2S N1 = N1S C C CLOCKWISE BOUNDARY TRAVERSAL C 20 INDX = IEND(N1) - 1 NEXT = IADJ(INDX) IF ( DET(X(NEXT),Y(NEXT),Z(NEXT),X(N1),Y(N1),Z(N1), . XP,YP,ZP) .LT. 0. ) GO TO 21 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 22 IF (X(NEXT)*Q(1) + Y(NEXT)*Q(2) + Z(NEXT)*Q(3) . .GE. 0.) GO TO 22 C C P, NEXT, N1, AND N2 ARE NEARLY COLLINEAR AND N1 IS THE C RIGHTMOST VISIBLE NODE C NF = N1 21 N2 = N1 N1 = NEXT GO TO 20 C C N1 IS THE LEFTMOST VISIBLE NODE C 22 NL = N1 C C NF AND NL HAVE BEEN FOUND C 23 I1 = NF I2 = NL I3 = 0 RETURN C C ALL POINTS ARE COLLINEAR C 24 I1 = 0 I2 = 0 I3 = 0 RETURN C C MORE THAN NRMAX RESTARTS. PRINT AN ERROR MESSAGE AND C TERMINATE PROCESSING. C 25 WRITE (LUN,100) NST, N0, N1, N2, N3, N4 100 FORMAT (1H1,37HERROR IN TRFIND -- POSSIBLE INFINITE , . 23HLOOP DUE TO F.P. ERROR.//1H , . 21HNST,N0,N1,N2,N3,N4 = ,5(I5,2H, ),I5//1H , . 40HTRY REORDERING THE NODES OR CHANGING NST) STOP END SUBROUTINE TRMESH (N,X,Y,Z, IADJ,IEND,IER) INTEGER N, IADJ(1), IEND(N), IER REAL X(N), Y(N), Z(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE CREATES A THIESSEN TRIANGULATION OF THE C CONVEX HULL OF N ARBITRARILY DISTRIBUTED POINTS (NODES) ON C THE UNIT SPHERE. IF THE NODES ARE NOT CONTAINED IN A C SINGLE HEMISPHERE, THEIR CONVEX HULL IS THE ENTIRE SPHERE, C AND THERE ARE NO BOUNDARY NODES. THE TRIANGULATION IS C OPTIMAL IN THE SENSE THAT IT IS AS NEARLY EQUIANGULAR AS C POSSIBLE. TRMESH IS PART OF AN INTERPOLATION PACKAGE C WHICH ALSO PROVIDES SUBROUTINES TO REORDER THE NODES, C TRANSFORM COORDINATES, ADD A NEW NODE, DELETE AN ARC, PLOT C THE MESH, AND PRINT THE DATA STRUCTURE. C UNLESS THE NODES ARE ALREADY ORDERED IN SOME REASONABLE C FASHION, THEY SHOULD BE REORDERED BY SUBROUTINE REORDR FOR C INCREASED EFFICIENCY BEFORE CALLING TRMESH. SPHERICAL C COORDINATES (LATITUDE AND LONGITUDE) MAY BE CONVERTED TO C CARTESIAN COORDINATES BY SUBROUTINE TRANS. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C N .GE. 3. C C X,Y,Z - N-VECTORS CONTAINING CARTESIAN C COORDINATES OF DISTINCT NODES. C (X(I),Y(I),Z(I)) DEFINES NODE C I. X(I)**2 + Y(I)**2 + Z(I)**2 C = 1 FOR ALL I. THE FIRST THREE C NODES MUST NOT BE COLLINEAR C (LIE ON A SINGLE GREAT CIRCLE). C C IADJ - VECTOR OF LENGTH .GE. 6*(N-1). C C IEND - VECTOR OF LENGTH .GE. N. C C N, X, Y, AND Z ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - IADJ - ADJACENCY LISTS OF NEIGHBORS IN C COUNTERCLOCKWISE ORDER. THE C LIST FOR NODE I+1 FOLLOWS THAT C FOR NODE I. THE VALUE 0 DENOTES C THE BOUNDARY (OR A PSEUDO-NODE C FROM WHICH ALL BOUNDARY NODES C ARE VISIBLE) AND IS ALWAYS THE C LAST NEIGHBOR OF A BOUNDARY C NODE. IADJ IS UNCHANGED IF IER C .NE. 0. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS (SETS OF C NEIGHBORS) IN IADJ. THE C NEIGHBORS OF NODE 1 BEGIN IN C IADJ(1). FOR K .GT. 1, THE C NEIGHBORS OF NODE K BEGIN IN C IADJ(IEND(K-1)+1) AND K HAS C IEND(K) - IEND(K-1) NEIGHBORS C INCLUDING (POSSIBLY) THE C BOUNDARY. IADJ(IEND(K)) .EQ. 0 C IFF NODE K IS ON THE BOUNDARY. C IEND IS UNCHANGED IF IER .NE. C 0. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF N .LT. 3. C IER = 2 IF THE FIRST THREE C NODES ARE COLLINEAR. C C MODULES REFERENCED BY TRMESH - ADNODE, TRFIND, INTADD, C BDYADD, COVSPH, SHIFTD, C INDEX, SWPTST, SWAP C C*********************************************************** C INTEGER NN, K, IERR LOGICAL LEFT C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C K = INDEX OF NODE TO BE ADDED, DO-LOOP INDEX C IERR = ERROR FLAG FOR CALL TO ADNODE - NOT CHECKED C LEFT = STATEMENT FUNCTION -- TRUE IFF (X0,Y0,Z0) IS IN C THE (CLOSED) LEFT HEMISPHERE DEFINED BY THE C PLANE CONTAINING (0,0,0), (X1,Y1,Z1), AND C (X2,Y2,Z2) WHERE LEFT IS DEFINED RELATIVE TO C AN OBSERVER AT (X1,Y1,Z1) FACING (X2,Y2,Z2) C LEFT (X1,Y1,Z1,X2,Y2,Z2,X0,Y0,Z0) = X0*(Y1*Z2-Y2*Z1) . - Y0*(X1*Z2-X2*Z1) + Z0*(X1*Y2-X2*Y1) .GE. 0. NN = N IER = 1 IF (NN .LT. 3) RETURN IER = 0 IF ( .NOT. LEFT (X(1),Y(1),Z(1),X(2),Y(2),Z(2), . X(3),Y(3),Z(3)) ) GO TO 2 IF ( .NOT. LEFT (X(2),Y(2),Z(2),X(1),Y(1),Z(1), . X(3),Y(3),Z(3)) ) GO TO 1 C C FIRST 3 NODES ARE COLLINEAR C IER = 2 RETURN C C FIRST TRIANGLE IS (1,2,3), I.E. 3 STRICTLY LEFT 1-2, C I.E. 3 LIES IN THE LEFT HEMISPHERE DEFINED BY ARC 1-2 C 1 IADJ(1) = 2 IADJ(2) = 3 IADJ(3) = 0 IADJ(4) = 3 IADJ(5) = 1 IADJ(6) = 0 IADJ(7) = 1 IADJ(8) = 2 IADJ(9) = 0 GO TO 3 C C FIRST TRIANGLE IS (3,2,1) = (2,1,3) = (1,3,2). INITIALIZE C IADJ. C 2 IADJ(1) = 3 IADJ(2) = 2 IADJ(3) = 0 IADJ(4) = 1 IADJ(5) = 3 IADJ(6) = 0 IADJ(7) = 2 IADJ(8) = 1 IADJ(9) = 0 C C INITIALIZE IEND C 3 IEND(1) = 3 IEND(2) = 6 IEND(3) = 9 IF (NN .EQ. 3) RETURN C C ADD NODES 4,...,N TO THE TRIANGULATION C DO 4 K = 4,NN CALL ADNODE (K,X,Y,Z, IADJ,IEND, IERR) 4 CONTINUE RETURN END SUBROUTINE TRPLOT (N,X,Y,Z,IADJ,IEND,ELAT,ELON, . ITITLE,NC,NUMBR, IER) INTEGER N, IADJ(1), IEND(N), ITITLE(1), NC, NUMBR, IER REAL X(N), Y(N), Z(N), ELAT, ELON C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE PLOTS THE ARCS OF A TRIANGULATION OF A C SET OF NODES ON THE UNIT SPHERE. VISIBLE NODES ARE PRO- C JECTED ONTO A PLANE WHICH CONTAINS THE ORIGIN AND HAS C NORMAL GIVEN BY A USER-SPECIFIED EYE-POSITION. PROJEC- C TIONS OF ADJACENT (VISIBLE) NODES ARE CONNECTED BY LINE C SEGMENTS AND A UNIT CIRCLE (THE PROJECTION OF THE SPHERE C BOUNDARY) IS DRAWN AROUND THE (PROJECTED) TRIANGULATION. C A SPECIAL SYMBOL IS PLOTTED IN THE POSITION OF A NODE C WHOSE NEIGHBORS ARE NOT ALL VISIBLE FROM THE SPECIFIED C EYE-POSITION. CARDS WITH C* IN THE FIRST TWO COLUMNS C MUST BE REPLACED WITH CALLS TO USER-SUPPLIED GRAPHICS C SUBROUTINES IN ORDER TO MAKE USE OF THIS ROUTINE. C C INPUT PARAMETERS - N - NUMBER OF NODES. N .GE. 3. C C X,Y,Z - N-VECTORS CONTAINING THE C CARTESIAN COORDINATES OF THE C NODES. C C IADJ,IEND - TRIANGULATION DATA STRUCTURE C (SEE SUBROUTINE TRMESH). C C ELAT,ELON - LATITUDE AND LONGITUDE OF THE C ENDPOINT OF THE UNIT NORMAL TO C THE PROJECTION PLANE. THE C VECTOR FROM THE ORIGIN (0,0,0) C TO THE EYE-POSITION IS TAKEN TO C BE A LARGE MULTIPLE OF THIS UNIT C NORMAL VECTOR. THUS ELAT AND C ELON SPECIFY THE POINT AT THE C CENTER OF THE PLOT. -PI/2 .LE. C ELAT .LE. PI/2 AND -PI .LE. ELON C .LE. PI. C C ITITLE - INTEGER ARRAY CONTAINING A LINE C OF TEXT TO BE CENTERED ABOVE THE C PLOT IF NC .GT. 0. ITITLE MUST C BE INITIALIZED WITH HOLLERITH C CONSTANTS OR READ WITH AN A-FOR- C MAT. ITS DIMENSION DEPENDS ON C NC AND THE NUMBER OF CHARACTERS C STORED IN A COMPUTER WORD. C C NC - NUMBER OF CHARACTERS IN ITITLE. C 0 .LE. NC .LE. 40. NO TITLE IS C DRAWN IF NC = 0. C C NUMBR - OPTION INDICATOR. IF NUMBR .NE. C 0, THE NODAL INDICES ARE PLOTTED C NEXT TO THE NODES. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE C ENCOUNTERED. C IER = 1 IF N OR NC IS OUT OF C RANGE. C C THE VALUES IN THE DATA STATEMENT BELOW MAY BE ALTERED C IN ORDER TO MODIFY VARIOUS PLOTTING OPTIONS. C C MODULES REFERENCED BY TRPLOT - CIRCLE, PRJCT C C INTRINSIC FUNCTIONS CALLED BY TRPLOT - COS, SIN, ATAN, C ABS, IABS C C*********************************************************** C INTEGER NPC, NN, ISW, N0, INDF, INDL, I, IERR, N1, NEW REAL ER, TOL, CX(481), CY(481), EX, EY, EZ, VX, VY, . VZ, X0, Y0, X1, Y1 LOGICAL ST0, BDRY DATA NPC, ER, TOL/ . 481, 50., .004/ C C LOCAL PARAMETERS - C C NPC = DIMENSION OF CX AND CY, AND NUMBER OF POINTS C (PLUS 1) USED TO PLOT THE UNIT CIRCLE -- C DETERMINES HOW CLOSELY THE CIRCLE IS APPROX- C IMATED. NOTE THAT IF NPC IS INCREASED, THE C DIMENSION STATEMENTS MUST ALSO BE INCREASED. C ER = RADIAL COORDINATE OF THE EYE-POSITION (DISTANCE C FROM THE ORIGIN) C TOL = MINIMUM ALLOWABLE ARC-LENGTH BETWEEN THE EYE- C POSITION AND (VX,VY,VZ) C NN = LOCAL COPY OF N C ISW = SWITCH USED BY SUBROUTINE PRJCT C N0 = NODE WHOSE INCIDENT ARCS ARE TO BE PLOTTED C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF N0 C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF N0 C I = IADJ INDEX IN THE RANGE INDF,...,INDL C IERR = ERROR FLAG FOR CALLS TO PRJCT -- NOT CHECKED C N1 = NEIGHBOR OF N0 C NEW = NEIGHBOR OF N0 AND CANDIDATE FOR NEXT VALUE OF C N0 C CX,CY = ARRAYS CONTAINING THE COORDINATES OF POINTS ON C THE UNIT CIRCLE C EX,EY,EZ = CARTESIAN COORDINATES OF THE EYE-POSITION C VX,VY,VZ = CARTESIAN COORDINATES OF A POINT WHOSE PROJECT- C ION DETERMINES THE POSITIVE Y-AXIS FOR THE C PLOT C X0,Y0 = PROJECTED COORDINATES OF N0 C X1,Y1 = PROJECTED COORDINATES OF N1 C ST0 = SWITCH USED TO ALTERNATE DIRECTION OF PEN C MOVEMENT C BDRY = SWITCH DETERMINING WHETHER ALL NEIGHBORS OF N0 C ARE VISIBLE C NN = N C C CHECK FOR INVALID PARAMETERS C IER = 1 IF (NN .LT. 3 .OR. NC .LT. 0 .OR. NC .GT. 40) RETURN IER = 0 C C COMMANDS WHICH PERFORM THE FOLLOWING FUNCTIONS SHOULD BE C INSERTED HERE -- C C* INITIALIZE THE PLOTTING ENVIRONMENT IF NECESSARY, C* SET DIMENSIONS OF THE PLOTTER SPACE, C* ESTABLISH A LINEAR MAPPING FROM THE DATA SPACE (-1,1) X C* (-1,1) TO THE PLOTTER SPACE, AND C* PLOT THE TITLE IF NC .NE. 0. C C PLOT THE UNIT CIRCLE C CALL CIRCLE(NPC, CX,CY) C* CALL CURVE(CX,CY,NPC) C C SET PARAMETERS FOR PRJCT. THE Y-AXIS WILL BE THE PRO- C JECTION OF THE NORTH POLE UNLESS THAT IS TOO CLOSE TO C THE EYE-POSITION. C EX = ER*COS(ELON)*COS(ELAT) EY = ER*SIN(ELON)*COS(ELAT) EZ = ER*SIN(ELAT) VX = 0. VY = 0. VZ = 1. IF (2.*ATAN(1.) - ABS(ELAT) .GE. TOL) GO TO 1 VX = 1. IF (EZ .GT. 0.) VX = -1. VZ = 0. 1 ISW = 1 C C INITIALIZE FOR LOOP ON NODES. EACH VISIBLE NODE N0 IS C CONNECTED TO ITS (VISIBLE) NEIGHBORS WHICH HAVE LARGER C INDICES. N0 IS THEN MARKED BY MAKING THE CORRESPONDING C IEND ENTRY NEGATIVE, AND THE SEARCH FOR THE NEXT C UNMARKED NODE BEGINS WITH THE NEIGHBORS OF N0. C N0 = 1 INDF = 1 C C TOP OF LOOP -- SET INDL AND TEST FOR N0 VISIBLE C 2 INDL = IEND(N0) IF (EX*X(N0) + EY*Y(N0) + EZ*Z(N0) .GE. 0.) GO TO 3 IEND(N0) = -INDL I = INDF GO TO 8 C C INITIALIZE ST0 AND BDRY, COMPUTE X0,Y0, AND INITIALIZE I C FOR NEIGHBOR LOOP C 3 ST0 = .TRUE. BDRY = .FALSE. CALL PRJCT(X(N0),Y(N0),Z(N0),0.,0.,0.,EX,EY,EZ, . VX,VY,VZ, ISW, X0,Y0,IERR) I = INDL IF (IADJ(I) .EQ. 0) I = I-1 C C LOOP ON NEIGHBORS OF N0 IN REVERSE ORDER C 4 N1 = IADJ(I) IF (EX*X(N1) + EY*Y(N1) + EZ*Z(N1) .GE. 0.) GO TO 5 C C N1 NOT VISIBLE -- SET BDRY AND MARK N1 C BDRY = .TRUE. IEND(N1) = -IABS(IEND(N1)) GO TO 6 C C N1 VISIBLE C 5 IF (N1 .LT. N0) GO TO 6 C C N1 IS VISIBLE AND HAS LARGER INDEX THAN N0 -- COMPUTE ITS C PROJECTION C CALL PRJCT(X(N1),Y(N1),Z(N1),0.,0.,0.,EX,EY,EZ, . VX,VY,VZ, ISW, X1,Y1,IERR) C C CONNECT N0 AND N1 -- THE DIRECTION OF PEN MOVEMENT ALTER- C NATES BETWEEN AWAY FROM N0 AND TOWARD N0 FOR REDUCED C PEN-UP TIME C C* IF (ST0) CALL LINE(X0,Y0,X1,Y1) C* IF (.NOT. ST0) CALL LINE(X1,Y1,X0,Y0) ST0 = .NOT. ST0 C C TEST FOR TERMINATION OF NEIGHBOR LOOP C 6 IF (I .EQ. INDF) GO TO 7 I = I - 1 GO TO 4 C C MARK N0 AS HAVING BEEN PROCESSED, PLOT A SYMBOL IF ANY OF C ITS NEIGHBORS ARE NOT VISIBLE, AND PLOT THE INDEX IF C NUMBR .NE. 0. C 7 IEND(N0) = -INDL C* IF (BDRY) CALL POINT(X0,Y0) C* IF (NUMBR .NE. 0) CALL RLINT(N0,X0,Y0) C C SEARCH THE NEIGHBORS OF N0 FOR AN UNMARKED NODE STARTING C WITH IADJ(I) = IADJ(INDF) C 8 NEW = IADJ(I) IF (NEW .EQ. 0) GO TO 10 IF (IEND(NEW) .LT. 0) GO TO 9 N0 = NEW INDF = IABS(IEND(N0-1)) + 1 GO TO 2 C C TEST FOR TERMINATION C 9 IF (I .EQ. INDL) GO TO 10 I = I + 1 GO TO 8 C C ALL NEIGHBORS OF N0 ARE MARKED. SEARCH IEND FOR AN C UNMARKED NODE C 10 DO 11 N0 = 2,NN IF (IEND(N0) .LT. 0) GO TO 11 INDF = -IEND(N0-1) + 1 GO TO 2 11 CONTINUE C C ALL NODES HAVE BEEN MARKED -- RESTORE IEND C DO 12 N0 = 1,NN 12 IEND(N0) = -IEND(N0) C C TERMINATE PLOTTING -- MOVE TO A NEW FRAME C C* CALL FRAME RETURN END SUBROUTINE TRPRNT (N,LUNIT,X,Y,Z,IADJ,IEND,IFLAG) INTEGER N, LUNIT, IADJ(1), IEND(N), IFLAG REAL X(N), Y(N), Z(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A TRIANGULATION OF A SET OF POINTS ON THE UNIT C SPHERE, THIS ROUTINE PRINTS THE ADJACENCY LISTS AND, C OPTIONALLY, THE NODAL COORDINATES (EITHER LATITUDE AND C LONGITUDE OR CARTESIAN COORDINATES). THE NUMBERS OF C BOUNDARY NODES, TRIANGLES, AND ARCS ARE ALSO PRINTED. C C INPUT PARAMETERS - N - NUMBER OF NODES IN THE MESH. C 3 .LE. N .LE. 9999. C C LUNIT - LOGICAL UNIT FOR OUTPUT. 1 C .LE. LUNIT .LE. 99. OUTPUT IS C PRINTED ON UNIT 6 IF LUNIT IS C OUT OF RANGE. C C X,Y,Z - VECTORS OF COORDINATES OF THE C NODES IN THE MESH. THESE MAY C BE DUMMY PARAMETERS IF IFLAG IS C NOT 0 OR 1, AND Z IS NOT USED C IF IFLAG .NE. 0. C C IADJ - SET OF ADJACENCY LISTS OF NODES C IN THE MESH. C C IEND - POINTERS TO THE ENDS OF C ADJACENCY LISTS IN IADJ FOR C EACH NODE IN THE MESH. C C IFLAG - OPTION INDICATOR C IFLAG = 0 IF X, Y, AND Z ARE TO C BE PRINTED (TO 6 DEC- C IMAL PLACES). C IFLAG = 1 IF ONLY X AND Y (AS- C SUMED TO CONTAIN LON- C GITUDE AND LATITUDE, C RESPECTIVELY) ARE TO C BE PRINTED. C C IADJ AND IEND MAY BE CREATED BY TRMESH. C C NONE OF THE PARAMETERS ARE ALTERED BY THIS ROUTINE. C C MODULES REFERENCED BY TRPRNT - NONE C C*********************************************************** C INTEGER NN, NMAX, LUN, NODE, INDF, INDL, NL, NLMAX, . INC, I, NB, NT, NA DATA NMAX/9999/, NLMAX/60/ C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C NMAX = UPPER BOUND ON N C LUN = LOCAL COPY OF LUNIT C NODE = INDEX OF A NODE C INDF,INDL = IADJ INDICES OF THE FIRST AND LAST NEIGHBORS C OF NODE C NL = NUMBER OF LINES PRINTED ON A PAGE C NLMAX = MAXIMUM NUMBER OF PRINT LINES PER PAGE EXCEPT C FOR THE LAST PAGE WHICH MAY HAVE 2 ADDITION- C AL LINES C INC = INCREMENT FOR NL C I = IADJ INDEX FOR IMPLIED DO-LOOP C NB = NUMBER OF BOUNDARY NODES C NT = NUMBER OF TRIANGLES C NA = NUMBER OF ARCS (UNDIRECTED EDGES) C NN = N LUN = LUNIT IF (LUN .LT. 1 .OR. LUN .GT. 99) LUN = 6 C C PRINT HEADING AND INITIALIZE NL C WRITE (LUN,100) NN IF (NN .LT. 3 .OR. NN .GT. NMAX) GO TO 7 NL = 6 IF (IFLAG .EQ. 0) GO TO 4 IF (IFLAG .EQ. 1) GO TO 2 C C PRINT IADJ ONLY C WRITE (LUN,101) NB = 0 INDF = 1 DO 1 NODE = 1,NN INDL = IEND(NODE) IF (IADJ(INDL) .EQ. 0) NB = NB + 1 INC = (INDL-INDF)/14 + 2 NL = NL + INC IF (NL .GT. NLMAX) WRITE (LUN,108) IF (NL .GT. NLMAX) NL = INC WRITE (LUN,104) NODE, (IADJ(I), I = INDF,INDL) IF (INDL-INDF .NE. 13) WRITE (LUN,107) 1 INDF = INDL + 1 GO TO 6 C C PRINT X (LONGITUDE), Y (LATITUDE), AND IADJ C 2 WRITE (LUN,102) NB = 0 INDF = 1 DO 3 NODE = 1,NN INDL = IEND(NODE) IF (IADJ(INDL) .EQ. 0) NB = NB + 1 INC = (INDL-INDF)/8 + 2 NL = NL + INC IF (NL .GT. NLMAX) WRITE (LUN,108) IF (NL .GT. NLMAX) NL = INC WRITE (LUN,105) NODE, X(NODE), Y(NODE), . (IADJ(I), I = INDF,INDL) IF (INDL-INDF .NE. 7) WRITE (LUN,107) 3 INDF = INDL + 1 GO TO 6 C C PRINT X, Y, Z, AND IADJ C 4 WRITE (LUN,103) NB = 0 INDF = 1 DO 5 NODE = 1,NN INDL = IEND(NODE) IF (IADJ(INDL) .EQ. 0) NB = NB + 1 INC = (INDL-INDF)/5 + 2 NL = NL + INC IF (NL .GT. NLMAX) WRITE (LUN,108) IF (NL .GT. NLMAX) NL = INC WRITE (LUN,106) NODE, X(NODE), Y(NODE), Z(NODE), . (IADJ(I), I = INDF,INDL) IF (INDL-INDF .NE. 4) WRITE (LUN,107) 5 INDF = INDL + 1 C C PRINT NB, NA, AND NT C 6 NT = 2*NN - NB - 2 IF (NB .EQ. 0) NT = NT - 2 NA = NT + NN - 1 IF (NB .EQ. 0) NA = NA - 1 WRITE (LUN,109) NB, NA, NT RETURN C C N IS OUT OF RANGE C 7 WRITE (LUN,110) RETURN C C PRINT FORMATS C 100 FORMAT (1H1,26X,23HADJACENCY SETS, N = ,I5//) 101 FORMAT (1H ,4HNODE,32X,17HNEIGHBORS OF NODE//) 102 FORMAT (1H ,4HNODE,4X,9HLONGITUDE,7X,8HLATITUDE, . 19X,17HNEIGHBORS OF NODE//) 103 FORMAT (1H ,4HNODE,5X,7HX(NODE),8X,7HY(NODE),8X, . 7HZ(NODE),12X,17HNEIGHBORS OF NODE//) 104 FORMAT (1H ,I4,5X,14I5/(1H ,9X,14I5)) 105 FORMAT (1H ,I4,2E15.6,5X,8I5/(1H ,39X,8I5)) 106 FORMAT (1H ,I4,3E15.6,5X,5I5/(1H ,54X,5I5)) 107 FORMAT (1H ) 108 FORMAT (1H1) 109 FORMAT (/1H ,5HNB = ,I4,15H BOUNDARY NODES,10X, . 5HNA = ,I5,5H ARCS,10X,5HNT = ,I5, . 10H TRIANGLES) 110 FORMAT (1H ,10X,25H*** N IS OUT OF RANGE ***) END SUBROUTINE UNIF (N,X,Y,Z,W,IADJ,IEND,NROW,NI,NJ,PLAT, . PLON,IFLAG, GRAD, WW,IER) INTEGER N, IADJ(1), IEND(N), NROW, NI, NJ, IFLAG, IER REAL X(N), Y(N), Z(N), W(N), PLAT(NI), PLON(NJ), . GRAD(3,N), WW(NROW,NJ) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A THIESSEN TRIANGULATION OF A SET OF NODES ON THE C UNIT SPHERE AND CORRESPONDING DATA VALUES, THIS ROUTINE C INTERPOLATES THE DATA VALUES TO A UNIFORM GRID FOR SUCH C APPLICATIONS AS CONTOURING. THE INTERPOLANT IS ONCE CON- C TINUOUSLY DIFFERENTIABLE. EXTRAPOLATION IS PERFORMED AT C GRID POINTS EXTERIOR TO THE TRIANGULATION WHEN THE NODES C DO NOT COVER THE ENTIRE SPHERE. C C INPUT PARAMETERS - N - NUMBER OF NODES. N .GE. 3 AND C N .GE. 7 IF IFLAG .NE. 1. C C X,Y,Z - VECTORS CONTAINING CARTESIAN COORD- C INATES OF THE NODES. C C W - VECTOR CONTAINING DATA VALUES. C W(I) IS ASSOCIATED WITH (X(I), C Y(I),Z(I)). C C IADJ,IEND - TRIANGULATION DATA STRUCTURE - MAY C BE CREATED BY TRMESH. C C NROW - NUMBER OF ROWS IN THE DIMENSION C STATEMENT OF WW. C C NI,NJ - NUMBER OF ROWS AND COLUMNS IN THE C UNIFORM GRID. 1 .LE. NI .LE. NROW, C 1 .LE. NJ. C C PLAT,PLON - VECTORS OF LENGTH NI AND NJ, RE- C SPECTIVELY, CONTAINING THE LATI- C TUDES AND LONGITUDES OF THE GRID C LINES. C C IFLAG - OPTION INDICATOR C IFLAG = 0 IF GRADIENT ESTIMATES AT C THE VERTICES OF A TRIAN- C GLE ARE TO BE RECOMPUTED C FOR EACH GRID POINT IN C THE TRIANGLE AND NOT C SAVED. C IFLAG = 1 IF GRADIENT ESTIMATES ARE C INPUT IN GRAD. C IFLAG = 2 IF GRADIENT ESTIMATES ARE C TO BE COMPUTED ONCE FOR C EACH NODE (BY GRADL) AND C SAVED IN GRAD. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C GRAD - NOT USED IF IFLAG = 0, 3 BY N ARRAY C WHOSE COLUMNS CONTAIN THE X, Y, AND C Z COMPONENTS (IN THAT ORDER) OF THE C GRADIENTS AT THE NODES IF IFLAG = C 1, OR ARRAY OF SUFFICIENT SIZE IF C IFLAG = 2. C C GRADIENT ESTIMATES MAY BE COMPUTED BY GRADL OR GRADG IF C IFLAG = 1. C C WW - NROW BY NCOL ARRAY WITH NROW .GE. C NI AND NCOL .GE. NJ. C C OUTPUT PARAMETERS - GRAD - ARRAY CONTAINING ESTIMATED C GRADIENTS AS DEFINED ABOVE IF C IFLAG = 2 AND IER .GE. 0. GRAD C IS NOT ALTERED IF IFLAG .NE. 2. C C WW - INTERPOLATED VALUES AT THE GRID C POINTS IF IER .GE. 0. WW(I,J) C = F(PLAT(I),PLON(J)) FOR I = C 1,...,NI AND J = 1,...,NJ. C C IER - ERROR INDICATOR C IER = K IF IF NO ERRORS WERE C ENCOUNTERED AND K GRID C POINTS REQUIRED EX- C TRAPOLATION FOR K .GE. C 0. C IER = -1 IF N, NI, NJ, OR IFLAG C IS OUT OF RANGE. C IER = -2 IF THE NODES ARE COL- C LINEAR. C IER = -3 IF EXTRAPOLATION FAIL- C ED DUE TO THE UNIFORM C GRID EXTENDING TOO FAR C BEYOND THE TRIANGULA- C TION BOUNDARY. C C MODULES REFERENCED BY UNIF - INTRC1, TRFIND, WVAL, ARCINT, C ARCLEN, C (AND OPTIONALLY) GRADL, GETNP, CONSTR, APLYR, C SETUP, GIVENS, ROTATE, APLYRT C C*********************************************************** C INTEGER NST, IST, NN, NX, NY, IFL, I, J, IERR, NEX DATA NST/1/ C C LOCAL PARAMETERS - C C NST = INITIAL VALUE FOR IST C IST = PARAMETER FOR INTRC1 C NN = LOCAL COPY OF N C NX,NY = LOCAL COPIES OF NI AND NJ C IFL = LOCAL COPY OF IFLAG C I,J = DO-LOOP INDICES C IERR = ERROR FLAG FOR CALLS TO GRADL AND INTRC1 C NEX = NUMBER OF GRID POINTS EXTERIOR TO THE TRIANGULA- C TION BOUNDARY (NUMBER OF EXTRAPOLATED VALUES) C NN = N NX = NI NY = NJ IFL = IFLAG IF (NX .LT. 1 .OR. NX .GT. NROW .OR. NY .LT. 1 . .OR. IFL .LT. 0 .OR. IFL .GT. 2) GO TO 4 IST = NST IF (IFL .NE. 2) GO TO 2 C C COMPUTE GRADIENT ESTIMATES AT THE NODES C DO 1 I = 1,NN CALL GRADL(NN,I,X,Y,Z,W,IADJ,IEND, GRAD(1,I),IERR) IF (IERR .LT. 0) GO TO 5 1 CONTINUE IFL = 1 C C COMPUTE UNIFORM GRID POINTS AND INTERPOLATED VALUES C 2 NEX = 0 DO 3 J = 1,NY DO 3 I = 1,NX CALL INTRC1(NN,PLAT(I),PLON(J),X,Y,Z,W,IADJ,IEND, . IFL,GRAD, IST, WW(I,J),IERR) IF (IERR .LT. 0) GO TO 5 3 NEX = NEX + IERR IER = NEX RETURN C C NI, NJ, OR IFLAG IS OUT OF RANGE C 4 IER = -1 RETURN C C ERROR IN GRADL OR INTRC1 C 5 IER = IERR RETURN END SUBROUTINE WVAL (B1,B2,B3,V1,V2,V3,W1,W2,W3,G1,G2,G3, . IFLAG, PW,PG) INTEGER IFLAG REAL B1, B2, B3, V1(3), V2(3), V3(3), W1, W2, W3, . G1(3), G2(3), G3(3), PW, PG(3) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN DATA VALUES AND GRADIENTS AT THE THREE VERTICES OF C A SPHERICAL TRIANGLE CONTAINING A POINT P, THIS ROUTINE C COMPUTES THE VALUE AND, OPTIONALLY, THE GRADIENT OF F AT P C WHERE F INTERPOLATES THE VERTEX DATA. ALONG THE TRIANGLE C EDGES, THE INTERPOLATORY FUNCTION F IS THE HERMITE CUBIC C (WITH RESPECT TO ARC-LENGTH) INTERPOLANT OF THE VALUES AND C TANGENTIAL GRADIENT COMPONENTS AT THE ENDPOINTS, AND THE C DERIVATIVE NORMAL TO THE ARC VARIES LINEARLY WITH RESPECT C TO ARC-LENGTH BETWEEN THE NORMAL GRADIENT COMPONENTS AT C THE ENDPOINTS. THUS THE METHOD YIELDS C-1 CONTINUITY WHEN C USED TO INTERPOLATE OVER A TRIANGULATION. THE INTERPOLANT C USES A FIRST-ORDER BLENDING METHOD ON THE UNDERLYING C PLANAR TRIANGLE. C C INPUT PARAMETERS - B1,B2,B3 - BARYCENTRIC COORDINATES OF C PP WITH RESPECT TO THE C (PLANAR) UNDERLYING TRIANGLE C (V1,V2,V3) WHERE PP IS THE C CENTRAL PROJECTION OF P ONTO C THIS TRIANGLE. C C V1,V2,V3 - CARTESIAN COORDINATES OF THE C VERTICES OF A SPHERICAL TRI- C ANGLE CONTAINING P. V3 LEFT C V1->V2. C C W1,W2,W3 - DATA VALUES ASSOCIATED WITH C THE VERTICES. C C G1,G2,G3 - GRADIENTS ASSOCIATED WITH C THE VERTICES. GI IS ORTHOG- C ONAL TO VI FOR I = 1,2,3. C C IFLAG - OPTION INDICATOR C IFLAG = 0 IF ONLY PW IS TO C BE COMPUTED. C IFLAG = 1 IF BOTH PW AND PG C ARE DESIRED. C C PG - VECTOR OF LENGTH 3 IF IFLAG C = 1, NOT USED OTHERWISE. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETERS - PW - INTERPOLATED VALUE AT P. C C PG - INTERPOLATED GRADIENT AT P C (ORTHOGONAL TO P) IF IFLAG C = 1. C C EACH VECTOR V ABOVE CONTAINS X-, Y-, AND Z-COMPONENTS IN C V(1), V(2), AND V(3), RESPECTIVELY. C C MODULES REFERENCED BY WVAL - ARCINT, ARCLEN C C INTRINSIC FUNCTION CALLED BY WVAL - SQRT C C*********************************************************** C INTEGER I REAL C1, C2, C3, SUM, U1(3), U2(3), U3(3), U1N, . U2N, U3N, Q1(3), Q2(3), Q3(3), VAL, W, G(3), . DUM C C LOCAL PARAMETERS - C C I = DO-LOOP INDEX C C1,C2,C3 = COEFFICIENTS (WEIGHT FUNCTIONS) OF PARTIAL C INTERPOLANTS. C1 = 1 ON THE EDGE OPPOSITE C V1 AND C1 = 0 ON THE OTHER EDGES. SIMI- C LARLY FOR C2 AND C3. C1+C2+C3 = 1. C SUM = QUANTITY USED TO NORMALIZE C1, C2, AND C3 C U1,U2,U3 = POINTS ON THE BOUNDARY OF THE PLANAR TRIAN- C GLE AND LYING ON THE LINES CONTAINING PP C AND THE VERTICES. U1 IS OPPOSITE V1, ETC. C U1N,U2N,U3N = QUANTITIES USED TO COMPUTE Q1, Q2, AND Q3 C (MAGNITUDES OF U1, U2, AND U3) C Q1,Q2,Q3 = CENTRAL PROJECTIONS OF U1, U2, AND U3 ONTO C THE SPHERE AND THUS LYING ON AN ARC OF THE C SPHERICAL TRIANGLE C VAL = LOCAL VARIABLE USED TO ACCUMULATE THE CON- C TRIBUTIONS TO PW C W,G = VALUE AND GRADIENT AT Q1, Q2, OR Q3 OBTAINED C BY INTERPOLATION ALONG ONE OF THE ARCS OF C THE SPHERICAL TRIANGLE C DUM = DUMMY VARIABLE FOR THE CALL TO ARCINT C C C COMPUTE WEIGHT FUNCTIONS C1, C2, AND C3 C C1 = B2*B3 C2 = B3*B1 C3 = B1*B2 SUM = C1 + C2 + C3 IF (SUM .GT. 0.) GO TO 2 C C P COINCIDES WITH A VERTEX C PW = B1*W1 + B2*W2 + B3*W3 IF (IFLAG .NE. 1) RETURN DO 1 I = 1,3 1 PG(I) = B1*G1(I) + B2*G2(I) + B3*G3(I) RETURN C C NORMALIZE C1, C2, AND C3 C 2 C1 = C1/SUM C2 = C2/SUM C3 = C3/SUM C C COMPUTE (U1,U2,U3) AND (U1N,U2N,U3N) C U1N = 0. U2N = 0. U3N = 0. DO 3 I = 1,3 U1(I) = (B2*V2(I) + B3*V3(I))/(B2+B3) U2(I) = (B3*V3(I) + B1*V1(I))/(B3+B1) U3(I) = (B1*V1(I) + B2*V2(I))/(B1+B2) U1N = U1N + U1(I)*U1(I) U2N = U2N + U2(I)*U2(I) 3 U3N = U3N + U3(I)*U3(I) C C COMPUTE Q1, Q2, AND Q3 C U1N = SQRT(U1N) U2N = SQRT(U2N) U3N = SQRT(U3N) DO 4 I = 1,3 Q1(I) = U1(I)/U1N Q2(I) = U2(I)/U2N 4 Q3(I) = U3(I)/U3N C C COMPUTE INTERPOLATED VALUE (VAL) AT P BY LOOPING ON C TRIANGLE SIDES C VAL = 0. C C CONTRIBUTION FROM SIDE OPPOSITE V1 -- C C COMPUTE VALUE AND GRADIENT AT Q1 BY INTERPOLATING C BETWEEN V2 AND V3 C CALL ARCINT (Q1,V2,V3,W2,W3,G2,G3, W,G,DUM) C C ADD IN THE CONTRIBUTION C VAL = VAL + C1*( W + B1*B1*(3.-2.*B1)*(W1-W) + . B1*(1.-B1)* . (B1*(G1(1)*U1(1)+G1(2)*U1(2)+G1(3)*U1(3)) . + (1.-B1)* . (G(1)*V1(1)+G(2)*V1(2)+G(3)*V1(3))/U1N) ) C C CONTRIBUTION FROM SIDE OPPOSITE V2 -- C C COMPUTE VALUE AND GRADIENT AT Q2 BY INTERPOLATING C BETWEEN V3 AND V1 C CALL ARCINT (Q2,V3,V1,W3,W1,G3,G1, W,G,DUM) C C ADD IN THE CONTRIBUTION C VAL = VAL + C2*( W + B2*B2*(3.-2.*B2)*(W2-W) + . B2*(1.-B2)* . (B2*(G2(1)*U2(1)+G2(2)*U2(2)+G2(3)*U2(3)) . + (1.-B2)* . (G(1)*V2(1)+G(2)*V2(2)+G(3)*V2(3))/U2N) ) C C CONTRIBUTION FROM SIDE OPPOSITE V3 -- C C COMPUTE INTERPOLATED VALUE AND GRADIENT AT Q3 C BY INTERPOLATING BETWEEN V1 AND V2 C CALL ARCINT (Q3,V1,V2,W1,W2,G1,G2, W,G,DUM) C C ADD IN THE FINAL CONTRIBUTION C VAL = VAL + C3*( W + B3*B3*(3.-2.*B3)*(W3-W) + . B3*(1.-B3)* . (B3*(G3(1)*U3(1)+G3(2)*U3(2)+G3(3)*U3(3)) . + (1.-B3)* . (G(1)*V3(1)+G(2)*V3(2)+G(3)*V3(3))/U3N) ) PW = VAL IF (IFLAG .NE. 1) RETURN C C COMPUTE PG C*** THIS AIN'T YET BEEN CODED *** C RETURN END C C C THIS DRIVER TESTS SOFTWARE PACKAGE SPHPAK FOR TRIANGULA- C TION AND INTERPOLATION ON THE SURFACE OF A SPHERE. ALL C MODULES OTHER THAN TRPLOT AND TRPRNT ARE TESTED UNLESS C A TRIANGULATION ERROR IS ENCOUNTERED, IN WHICH CASE THE C PROGRAM TERMINATES IMMEDIATELY. THE INTERPOLATION MODULES C ARE TESTED ON DATA TAKEN FROM A CONSTANT FUNCTION FOR C WHICH THE METHODS ARE EXACT, AND MAXIMUM ERRORS ARE C PRINTED. BY DEFAULT, TESTS ARE PERFORMED ON A SIMPLE C DATA SET CONSISTING OF 13 NODES -- THE NORTH POLE ALONG C WITH SIX NODES UNIFORMLY DISTRIBUTED AROUND THE 30-DEGREE C PARALLEL OF LATITUDE AND SIX NODES UNIFORMLY DISTRIBUTED C AROUND THE 60-DEGREE PARALLEL. HOWEVER, BY ENABLING THE C READ STATEMENTS BELOW (C* IN THE FIRST TWO COLUMNS), TEST- C ING MAY BE PERFORMED ON AN ARBITRARY SET OF UP TO 100 C NODES. C C THE I/O UNITS MAY BE CHANGED BY ALTERING LIN (INPUT) AND C LOUT (OUTPUT) IN THE DATA STATEMENT BELOW. IF NMAX IS C INCREASED, DIMENSION STATEMENTS MUST BE CHANGED ACCORDING- C LY. C INTEGER IADJ(600), IEND(100), NODES(200) REAL X(100), Y(100), Z(100), W(100), GRAD(3,100), . G(3), PLAT(4), PLON(5), WW(4,5) DATA LIN/5/, LOUT/6/, NMAX/100/, N/13/, . NROW/4/, NI/4/, NJ/5/ DATA X(1), X(2), X(3), X(4), X(5), X(6), X(7) . / 0., 0., 60., 120., 180., 240., 300./, . X(8), X(9), X(10), X(11), X(12), X(13) . / 30., 90., 150., 210., 270., 330./, . Y(1), Y(2), Y(3), Y(4), Y(5), Y(6), Y(7) . / 90., 30., 30., 30., 30., 30., 30./, . Y(8), Y(9), Y(10), Y(11), Y(12), Y(13) . / 60., 60., 60., 60., 60., 60./ DATA PLAT(1), PLAT(2), PLAT(3), PLAT(4) . / 0., 25., 50., 75./ DATA PLON(1), PLON(2), PLON(3), PLON(4), PLON(5) . / 0., 72., 144., 216., 288./ C C CONVERT DEGREES TO RADIANS C SC = ATAN(1.)/45. DO 1 K = 1,N X(K) = SC*X(K) 1 Y(K) = SC*Y(K) DO 2 I = 1,NI 2 PLAT(I) = SC*PLAT(I) DO 3 J = 1,NJ 3 PLON(J) = SC*PLON(J) C C INPUT THE NUMBER OF NODES AND THEIR COORDINATES (X AND Y C ARE LONGITUDE AND LATITUDE, RESPECTIVELY, IN RADIANS). C THE FIRST THREE NODES MUST NOT LIE ON A COMMON GREAT C CIRCLE. C C* READ (LIN,100) N IF (N .LT. 3 .OR. N .GT. NMAX) WRITE (LOUT,200) N IF (N .LT. 3 .OR. N .GT. NMAX) STOP C* READ (LIN,110) (X(I),Y(I), I = 1,N) 100 FORMAT (I4) 110 FORMAT (F11.8) C C PRINT A HEADING AND CONVERT SPHERICAL TO CARTESIAN C COORDINATES. C WRITE (LOUT,400) N CALL TRANS (N,Y,X, X,Y,Z) C C TEST TRMESH (REFER TO TRMTST) C CALL TRMESH (N,X,Y,Z, IADJ,IEND,IER) IF (IER .EQ. 2) WRITE (LOUT,210) IF (IER .EQ. 2) STOP CALL TRMTST (N,IADJ,IEND,LOUT, IER) IF (IER .GT. 0) STOP C C TEST GETNP BY ORDERING THE NODES ON DISTANCE FROM N0 AND C VERIFYING THE ORDERING. THE STUFF WITH KSUM IS C TOO CLEVER FOR WORDS. C N0 = N/2 NODES(1) = N0 KSUM = N0 DKM1 = -1.1 DO 4 K = 2,N CALL GETNP (X,Y,Z,IADJ,IEND,K, NODES, DK,IER) IF (IER .NE. 0 .OR. DK .LT. DKM1) GO TO 5 KSUM = KSUM + NODES(K) 4 DKM1 = DK IF (KSUM .EQ. (N*(N+1))/2) GO TO 6 5 WRITE (LOUT,220) STOP C C SET DATA VALUES AND GRADIENTS FOR INTERPOLATION TESTS C ON THE FUNCTION F(X,Y,Z) = 1. C 6 DO 7 K = 1,N W(K) = 1. GRAD(1,K) = 0. GRAD(2,K) = 0. 7 GRAD(3,K) = 0. C C TEST GRADL C DMAX = 0. DO 8 K = 1,N CALL GRADL (N,K,X,Y,Z,W,IADJ,IEND, G,IER) IF (IER .LT. 6) GO TO 9 8 DMAX = AMAX1(ABS(G(1)),ABS(G(2)),ABS(G(3)),DMAX) WRITE (LOUT,300) DMAX GO TO 10 C 9 WRITE (LOUT,310) IER STOP C C TEST GRADG C 10 TOL = 0. MAXIT = 1 CALL GRADG (N,X,Y,Z,W,IADJ,IEND,TOL, MAXIT,GRAD, IER) DMAX = 0. DO 11 K = 1,N 11 DMAX = AMAX1(ABS(GRAD(1,K)),ABS(GRAD(2,K)), . ABS(GRAD(3,K)),DMAX) WRITE (LOUT,320) DMAX C C TEST INTRC0 C DMAX = 0. IST = 1 NEX = 0 NFAIL = 0 DO 12 J = 1,NJ DO 12 I = 1,NI CALL INTRC0 (N,PLAT(I),PLON(J),X,Y,Z,W,IADJ,IEND, . IST, PW,IER) IF (IER .EQ. 1) NEX = NEX + 1 IF (IER .LT. 0) NFAIL = NFAIL + 1 IF (IER .LT. 0) GO TO 12 DMAX = AMAX1(DMAX,ABS(PW-1.)) 12 CONTINUE WRITE (LOUT,330) DMAX, NEX IF (NFAIL .GT. 0) WRITE (LOUT,410) NFAIL IF (NFAIL .GT. 0) GO TO 14 C C TEST UNIF (AND INTRC1) C CALL UNIF (N,X,Y,Z,W,IADJ,IEND,NROW,NI,NJ,PLAT, . PLON,1,GRAD, WW,IER) DMAX = 0. DO 13 J = 1,NJ DO 13 I = 1,NI 13 DMAX = AMAX1(ABS(WW(I,J)-1.),DMAX) WRITE (LOUT,340) DMAX, IER C C TEST REORDR (AND PERMUT) BY SORTING ON X-COMPONENTS AND C VERIFYING THE ORDERING. C 14 CALL REORDR (N,3, X,Y,Z,DUMMY, IEND) DO 15 K = 2,N IF (X(K) .LT. X(K-1)) GO TO 16 15 CONTINUE C C CREATE A NEW TRIANGULATION C CALL TRMESH (N,X,Y,Z, IADJ,IEND,IER) GO TO 17 16 WRITE (LOUT,230) STOP C C TEST EDGE, USING NODES AS WORK SPACE -- C FORCE AN EDGE BETWEEN NODE 1 AND A NODE WHICH C IS NOT ALREADY ADJACENT TO NODE 1 (IF IT C EXISTS). C 17 IF (IEND(1) .GE. N-1) WRITE (LOUT,420) IF (IEND(1) .GE. N-1) GO TO 21 N1 = 1 NP2 = N+2 DO 18 K = 2,N N2 = NP2 - K LWK = N CALL EDGE (N1,N2,X,Y,Z, LWK,NODES,IADJ,IEND, IER) IF (IER .NE. 0) GO TO 19 IF (LWK .GT. 0) GO TO 20 18 CONTINUE 19 WRITE (LOUT,240) IER STOP C C TEST FOR ALL ARCS (OTHER THAN N1-N2) OPTIMAL. C 20 L = 1 CALL ARCTST (N,X,Y,Z,IADJ,IEND, L, NODES,IER) IF (L .EQ. 1 .AND. NODES(1) .EQ. N1 .AND. . NODES(2) .EQ. N2) GO TO 21 IF (L .EQ. 0 .AND. LOPTST(N1,N2,X,Y,Z,IADJ,IEND) . .EQ. 0) GO TO 21 WRITE (LOUT,250) L STOP C C TEST DELETE BY REMOVING AS MANY BOUNDARY ARCS AS POSSIBLE C AND TESTING THE DATA STRUCTURE FOR ERRORS. C 21 NDL = 0 CALL BNODES (N,IADJ,IEND, NB,NA,NT,NODES) IF (NB .EQ. 0) WRITE (LOUT,430) IF (NB .EQ. 0) GO TO 25 N1 = NODES(NB) DO 22 K = 1,NB N2 = N1 N1 = NODES(K) CALL DELETE (N,N1,N2, IADJ,IEND, IER) IF (IER .EQ. 1 .OR. IER .EQ. 3) GO TO 23 22 IF (IER .EQ. 0) NDL = NDL + 1 IF (NDL .GT. 0) GO TO 24 WRITE (LOUT,430) GO TO 25 C 23 WRITE (LOUT,260) IER STOP C C CALL TRMTST C 24 CALL TRMTST (N,IADJ,IEND,0, IER) IF (IER .EQ. 0) GO TO 25 WRITE (LOUT,270) IER STOP C C SUCCESSFUL TEST C 25 WRITE (LOUT,440) STOP C C TRIANGULATION ERROR MESSAGE FORMATS C 200 FORMAT (//1H ,21HN OUT OF RANGE -- N =,I4) 210 FORMAT (1H ,23HALL NODES ARE COLLINEAR) 220 FORMAT (1H ,14HERROR IN GETNP) 230 FORMAT (1H ,15HERROR IN REORDR) 240 FORMAT (1H ,23HERROR IN EDGE -- IER = ,I1) 250 FORMAT (1H ,16HERROR IN EDGE --,I3, . 29H ARCS ARE NOT LOCALLY OPTIMAL) 260 FORMAT (1H ,30HERROR IN DELETE (INVALID FLAG , . 19HRETURNED) -- IER = ,I1) 270 FORMAT (1H ,31HERROR IN DELETE (TRIANGULATION , . 9HDESTROYED/1H ,25HERROR FLAG FROM TRMTST = , . I1) C C INTERPOLATION TEST MESSAGE FORMATS C 300 FORMAT (1H ,30HGRADL TEST -- MAXIMUM ERROR = , . E15.8//) 310 FORMAT (1H ,24HERROR IN GRADL -- IER = ,I1) 320 FORMAT (1H ,30HGRADG TEST -- MAXIMUM ERROR = , . E15.8//) 330 FORMAT (1H ,31HINTRC0 TEST -- MAXIMUM ERROR = , . E15.8/1H ,15X,17HEXTRAPOLATION AT ,I2, . 7H POINTS//) 340 FORMAT (1H ,30HUNIF (INTRC1) TEST -- MAXIMUM , . 8HERROR = ,E15.8/1H ,22X, . 26HEXTRAPOLATION OCCURRED AT ,I2,7H POINTS//) C C INFORMATIVE MESSAGE FORMATS C 400 FORMAT (1H1,17X,18HSPHPAK TEST -- N =,I4///) 410 FORMAT (1H ,38HEXTRAPOLATION IN INTRC0 TEST FAILED AT, . I3,10H POINTS --/1H , . 26HUNIF AND INTRC1 NOT TESTED//) 420 FORMAT (1H ,29HSUBROUTINE EDGE NOT TESTED --/ . 1H ,37H THE NODE WITH MINIMUM X-COORDINATE , . 30HIS ADJACENT TO ALL OTHER NODES//) 430 FORMAT (1H ,31HSUBROUTINE DELETE NOT TESTED --/ . 1H ,33H NO BOUNDARY ARC CAN BE REMOVED , . 7HWITHOUT/ . 1H ,30H DESTROYING THE TRIANGULATION//) 440 FORMAT (1H ,35HNO TRIANGULATION ERRORS ENCOUNTERED) END SUBROUTINE ARCTST (N,X,Y,Z,IADJ,IEND, LL, LIST,IER) INTEGER N, IADJ(1), IEND(N), LL, LIST(2,1), IER REAL X(N), Y(N), Z(N) LOGICAL SWPTST C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE DETERMINES THE NUMBER OF ARCS IN A TRIANGU- C LATION WHICH ARE NOT LOCALLY OPTIMAL AS DEFINED BY LOGICAL C FUNCTION SWPTST. A THIESSEN TRIANGULATION (SEE SUBROUTINE C TRMESH) IS CHARACTERIZED BY THE PROPERTY THAT ALL ARCS ARE C LOCALLY OPTIMAL. THE ALGORITHM CONSISTS OF APPLYING THE C SWAP TEST TO ALL INTERIOR ARCS. C C INPUT PARAMETERS - C C N - NUMBER OF NODES IN THE TRIANGULATION. N .GE. 3. C C X,Y,Z - COORDINATES OF THE NODES. C C IADJ,IEND - DATA STRUCTURE CONTAINING THE TRIANGULA- C TION. SEE SUBROUTINE TRMESH. C C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C LL - NUMBER OF COLUMNS RESERVED FOR LIST (SEE THE C OUTPUT VALUE OF LIST). LL .GE. 0. C C LIST - INTEGER ARRAY DIMENSIONED 2 BY LL (OR VECTOR C OF LENGTH .GE. 2*LL) IF LL .GT. 0. C C OUTPUT PARAMETERS - C C LL - NUMBER OF ARCS WHICH ARE NOT LOCALLY OPTIMAL C UNLESS IER = 1. C C LIST - COLUMNS CONTAIN THE ENDPOINT NODAL INDICES C (SMALLER INDEX IN THE FIRST ROW) OF THE FIRST C K NONOPTIMAL ARCS ENCOUNTERED, WHERE K IS THE C MINIMUM OF THE INPUT AND OUTPUT VALUES OF LL, C IF IER = 0 AND K .GT. 0. THE NUMBER OF INTE- C RIOR ARCS IS 3N-2NB-3 IF NB .NE. 0, AND 3N-6 C IF NB = 0 WHERE NB IS THE NUMBER OF BOUNDARY C NODES. BOUNDARY ARCS ARE OPTIMAL BY DEFINI- C TION. C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE ENCOUNTERED. C IER = 1 IF N .LT. 3 OR LL .LT. 0. C C MODULES REQUIRED BY ARCTST - SWPTST C C*********************************************************** C NM1 = N - 1 LMAX = LL C C TEST FOR ERRORS AND INITIALIZE FOR LOOP ON INTERIOR ARCS. C IF (NM1 .LT. 2 .OR. LMAX .LT. 0) GO TO 3 IER = 0 L = 0 INDF = 1 C C OUTER LOOP ON NODES IO1 C DO 2 IO1 = 1,NM1 INDL = IEND(IO1) IN2 = IADJ(INDL) C C INNER LOOP ON NEIGHBORS IO2 OF IO1 SUCH THAT IO2 .GT. IO1 C AND IO1-IO2 IS AN INTERIOR ARC -- (IO1,IO2,IN1) AND C (IO2,IO1,IN2) ARE TRIANGLES. C DO 1 INDX = INDF,INDL IO2 = IADJ(INDX) IN1 = IADJ(INDX+1) IF (INDX .EQ. INDL) IN1 = IADJ(INDF) IF (IO2 .LT. IO1 .OR. IN1 .EQ. 0 .OR. . IN2 .EQ. 0) GO TO 1 C C TEST FOR A SWAP. C IF (.NOT. SWPTST(IN1,IN2,IO1,IO2,X,Y,Z)) GO TO 1 L = L + 1 IF (L .GT. LMAX) GO TO 1 LIST(1,L) = IO1 LIST(2,L) = IO2 1 IN2 = IO2 2 INDF = INDL + 1 LL = L RETURN C C N OR LL OUT OF RANGE C 3 IER = 1 RETURN END INTEGER FUNCTION LOPTST (N1,N2,X,Y,Z,IADJ,IEND) INTEGER N1, N2, IADJ(1), IEND(1) REAL X(1), Y(1), Z(1) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C GIVEN A PAIR OF INDICES DEFINING A TRIANGULATION ARC, C THIS FUNCTION DETERMINES WHETHER OR NOT THE ARC IS LOCALLY C OPTIMAL AS DEFINED BY LOGICAL FUNCTION SWPTST. C C INPUT PARAMETERS - C C N1,N2 - X,Y,Z INDICES OF A PAIR OF ADJACENT NODES. C C X,Y,Z - NODAL COORDINATES. C C IADJ,IEND - DATA STRUCTURE CONTAINING THE TRIANGULA- C TION. SEE SUBROUTINE TRMESH. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - C C LOPTST = -2 IF N1 AND N2 ARE NOT ADJACENT, C = -1 IF N1-N2 IS AN INTERIOR ARC WHICH IS C NOT LOCALLY OPTIMAL, C = 0 IF N1-N2 SATISFIES THE NEUTRAL CASE (THE C VERTICES OF THE CORRESPONDING QUADRILAT- C ERAL LIE IN A COMMON PLANE), C = 1 IF N1-N2 IS A LOCALLY OPTIMAL INTERIOR C ARC, C = 2 IF N1-N2 IS A BOUNDARY ARC. C NOTE THAT N1-N2 IS LOCALLY OPTIMAL IFF LOPTST .GE. 0. C C MODULES REQUIRED BY LOPTST - NONE C C*********************************************************** C I1 = N1 I2 = N2 C C FIND THE INDEX OF I2 AS A NEIGHBOR OF I1 C INDF = 1 IF (I1 .GT. 1) INDF = IEND(I1-1) + 1 INDL = IEND(I1) 1 INDX = INDX - 1 IF (IADJ(INDX) .EQ. I2) GO TO 2 IF (INDX .NE. INDF) GO TO 1 C C N1 AND N2 ARE NOT ADJACENT C LOPTST = -2 RETURN C C DETERMINE I3 AND I4 SUCH THAT (I1,I2,I3) AND C (I2,I1,I4) ARE TRIANGLES. C 2 IF (INDX .NE. INDL) I3 = IADJ(INDX+1) IF (INDX .EQ. INDL) I3 = IADJ(INDF) IF (INDX .NE. INDF) I4 = IADJ(INDX-1) IF (INDX .EQ. INDF) I4 = IADJ(INDL) IF (I3 .NE. 0 .AND. I4 .NE. 0) GO TO 3 C C N1-N2 IS A BOUNDARY ARC C LOPTST = 2 RETURN C C COMPUTE COMPONENTS OF THE DIRECTED LINE SEGMENTS FROM C I4 TO THE OTHER THREE QUADRILATERAL VERTICES. C 3 DX1 = X(I1) - X(I4) DX2 = X(I2) - X(I4) DX3 = X(I3) - X(I4) DY1 = Y(I1) - Y(I4) DY2 = Y(I2) - Y(I4) DY3 = Y(I3) - Y(I4) DZ1 = Z(I1) - Z(I4) DZ2 = Z(I2) - Z(I4) DZ3 = Z(I3) - Z(I4) C C I1-I2 IS LOCALLY OPTIMAL IFF I3 LIES IN THE CLOSED HALF C SPACE DEFINED BY THE PLANE OF I1, I2, AND I4 AND WHICH C CONTAINS THE ORIGIN, I.E. IFF DET(I3-I4,I2-I4,I1-I4) C .LE. 0. C DET = DX3*(DY2*DZ1 - DY1*DZ2) . -DY3*(DX2*DZ1 - DX1*DZ2) . +DZ3*(DX2*DY1 - DX1*DY2) IF (DET .GT. 0.) GO TO 4 IF (DET .LT. 0.) GO TO 5 C C NEUTRAL CASE C LOPTST = 0 RETURN C C N1-N2 NOT LOCALLY OPTIMAL C 4 LOPTST = -1 RETURN C C N1-N2 LOCALLY OPTIMAL C 5 LOPTST = 1 RETURN END FUNCTION STORE (X) REAL X C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS FUNCTION FORCES ITS ARGMENT X TO BE STORED IN MAIN C MEMORY. THIS IS USEFUL FOR COMPUTING MACHINE DEPENDENT C PARAMETERS (SUCH AS THE MACHINE PRECISION) WHERE IT IS C NECESSARY TO AVOID COMPUTATION IN HIGH PRECISION REG- C ISTERS. C C INPUT PARAMETER - C C X - VALUE TO BE STORED. C C X IS NOT ALTERED BY THIS FUNCTION. C C OUTPUT PARAMETER - C C STORE - VALUE OF X AFTER IT HAS BEEN STORED AND C (POSSIBLY) TRUNCATED OR ROUNDED TO THE C MACHINE WORD LENGTH. C C MODULES REQUIRED BY STORE - NONE C C*********************************************************** C COMMON/STCOM/Y Y = X STORE = Y RETURN END SUBROUTINE TRMTST (N,IADJ,IEND,LUN, IER) INTEGER N, IADJ(1), IEND(N), LUN, IER C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE TESTS THE VALIDITY OF THE DATA STRUCTURE C REPRESENTING A THIESSEN TRIANGULATION CREATED BY SUBROU- C TINE TRMESH. THE FOLLOWING PROPERTIES ARE TESTED -- C 1) IEND(1) .GE. 3 AND IEND(K) .GE. IEND(K-1)+3 FOR K = C 2,...,N (EACH NODE HAS AT LEAST THREE NEIGHBORS). C 2) 0 .LE. IADJ(K) .LE. N FOR K = 1,...,IEND(N) (IADJ C ENTRIES ARE NODAL INDICES OR ZEROS REPRESENTING THE C BOUNDARY). C 3) THE TRIANGULATION PARAMETERS N, NB, NT, AND NA C (NUMBERS OF NODES, BOUNDARY NODES, TRIANGLES, AND C ARCS, RESPECTIVELY) SATISFY THE FOLLOWING RELATION- C SHIPS -- C NB .GE. 0, NB .NE. 1, AND NB .NE. 2, C NT = 2N-4 AND NA = 3N-6 IF NB = 0, AND C NT = 2N-NB-2 AND NA = 3N-NB-3 IF NB .GE. 3. C NO TEST IS MADE FOR THE PROPERTY THAT A TRIANGULATION C COVERS THE CONVEX HULL OF THE NODES OR THAT THE ARCS ARE C LOCALLY OPTIMAL. C C INPUT PARAMETERS - C C N - NUMBER OF NODES. N .GE. 3. C C IADJ,IEND - TRIANGULATION DATA STRUCTURE. SEE SUB- C ROUTINE TRMESH. C C LUN - LOGICAL UNIT NUMBER FOR PRINTING ERROR MES- C SAGES. IF LUN .LT. 1 OR LUN .GT. 99, NO MES- C SAGES ARE PRINTED. C C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - C C IER - ERROR INDICATOR C IER = 0 IF NO ERRORS WERE ENCOUNTERED. C IER = 1 IF N .LT. 3. C IER = 2 IF AN IEND OR IADJ ENTRY IS OUT OF C RANGE. C IER = 3 IF THE TRIANGULATION PARAMETERS ARE C INCONSISTENT. C C MODULES REQUIRED BY TRMTST - NONE C C*********************************************************** C LOGICAL RITE C C SET LOCAL VARIABLES, TEST FOR N OUT OF RANGE, AND C INITIALIZE A COUNT FOR NB. C NN = N RITE = LUN .GE. 1 .AND. LUN .LE. 99 IF (NN .LT. 3) GO TO 3 NB = 0 C C OUTER LOOP ON NODES N0 C INDF = 1 DO 2 N0 = 1,NN INDL = IEND(N0) IF (INDL .LT. INDF+2) GO TO 4 IF (IADJ(INDL) .EQ. 0) NB = NB + 1 C C INNER LOOP ON NEIGHBORS NAB OF N0 C DO 1 INDX = INDF,INDL NAB = IADJ(INDX) IF (NAB .LT. 0 .OR. NAB .GT. NN .OR. . (INDX .LT. INDL .AND. NAB .EQ. 0)) GO TO 5 1 CONTINUE 2 INDF = INDL + 1 C C CHECK PARAMETERS FOR CONSISTENCY -- NA = (L-NB)/2 AND C NT = (L-2*NB)/3 WHERE L IS THE LENGTH OF IADJ. C L = IEND(NN) IF ((NB .LT. 3 .AND. NB .NE. 0) .OR. . (NB .EQ. 0 .AND. L .NE. 6*NN-12) .OR. . (NB .NE. 0 .AND. L .NE. 6*NN-NB-6)) GO TO 6 C C NO ERRORS WERE ENCOUNTERED. C IER = 0 RETURN C C N IS OUT OF RANGE. C 3 IER = 1 IF (RITE) WRITE (LUN,100) N 100 FORMAT (1H ,32HTRMTST -- N OUT OF RANGE -- N = ,I5) RETURN C C IEND ENTRY OUT OF RANGE C 4 IER = 2 IF (RITE) WRITE (LUN,110) N0 110 FORMAT (1H ,15HTRMTST -- NODE ,I5, . 26H HAS LESS THAN 3 NEIGHBORS) RETURN C C IADJ ENTRY OUT OF RANGE C 5 IER = 2 IF (RITE) WRITE (LUN,120) INDX 120 FORMAT (1H ,33HTRMTST -- IADJ(K) IS NOT A VALID , . 14HINDEX FOR K = ,I5) RETURN C C INCONSISTENT TRIANGULATION PARAMETERS C 6 IER = 3 IF (.NOT. RITE) RETURN NA = (L-NB)/2 NT = (L-2*NB)/3 WRITE (LUN,130) NB, NT, NA 130 FORMAT (1H ,33HTRMTST -- INCONSISTENT PARAMETERS/ . 1H ,10X,I5,15H BOUNDARY NODES,3X,I5, . 10H TRIANGLES,3X,I5,5H ARCS) RETURN END C C C THIS DRIVER PROGRAM DEMONSTRATES USE OF THE PACKAGE C SPHPAK FOR INTERPOLATION OVER THE SURFACE OF A SPHERE. C A TRIANGULATION IS CREATED FROM A SET OF 25 NODES CON- C SISTING OF THE NORTH POLE AND 24 POINTS UNIFORMLY DIS- C TRIBUTED AROUND THE 60-DEGREE PARALLEL (WITH 15 DEGREES C LONGITUDINAL SEPARATION). THE TRIANGULATION DATA STRUC- C TURE IS PRINTED, AND A SET OF FUNCTION VALUES ON THE C NODES IS INTERPOLATED TO 100 POINTS UNIFORMLY DISTRIBUTED C OVER THE BANDED REGION CONTAINED BETWEEN 65 AND 85 DEGREES C LATITUDE. THE MAXIMUM AND MEAN DEVIATIONS FROM THE TRUE C FUNCTION VALUES AT THE INTERPOLATION POINTS ARE COMPUTED C AND PRINTED ON LOGICAL UNIT 6. THE FUNCTION CHOSEN C IS A CONSTANT FOR WHICH INTERPOLATION IS EXACT. C C FOR N NODES AND AN NI BY NJ UNIFORM GRID OF INTERPOLA- C TION POINTS, ARRAY DIMENSIONS ARE AS FOLLOWS -- C IADJ(6N-6), IEND(N), RLAT(N), RLON(N), X(N), Y(N), Z(N), C W(N), GRAD(3,N), PLAT(NI), PLON(NJ), AND WW(NROW,NJ) C FOR NROW .GE. NI. C INTEGER IADJ(144), IEND(25) REAL RLAT(25), RLON(25) REAL X(25), Y(25), Z(25), W(25), GRAD(3,25), . PLAT(5), PLON(20), WW(5,20) DATA N/25/, NROW/5/, NI/5/, NJ/20/ DATA LUNIT/6/ DATA OLAT, DLAT/65., 20./, OLON, DLON/0., 360./ TFUN(XF,YF,ZF) = 1. C C GENERATE THE SET OF NODES IN TERMS OF LATITUDINAL AND C LONGITUDINAL COORDINATES. THE SEPARATION DEL IS 15 C DEGREES. C RLAT(1) = 90. RLON(1) = 0. DEL = 360./FLOAT(N-1) DO 1 K = 2,N RLAT(K) = 60. 1 RLON(K) = FLOAT(K-2)*DEL C C SET X AND Y TO THE VALUES OF RLON AND RLAT, RESPECTIVELY, C IN RADIANS. (RLON AND RLAT ARE SAVED FOR PRINTING BY C TRPRNT). C SC = ATAN(1.)/45. DO 2 K = 1,N X(K) = SC*RLON(K) 2 Y(K) = SC*RLAT(K) C C TRANSFORM SPHERICAL COORDINATES X AND Y TO CARTESIAN C COORDINATES (X,Y,Z) ON THE UNIT SPHERE (X**2 + Y**2 + C Z**2 = 1). C CALL TRANS (N,Y,X, X,Y,Z) C C GENERATE DATA VALUES AT THE NODES FROM TFUN C DO 3 K = 1,N 3 W(K) = TFUN (X(K),Y(K),Z(K)) C C CREATE THE TRIANGULATION AND TEST THE ERROR FLAG C CALL TRMESH (N,X,Y,Z, IADJ,IEND,IER) IF (IER .NE. 0) WRITE (LUNIT,100) IER IF (IER .NE. 0) STOP C C PRINT THE SPHERICAL COORDINATES AND ADJACENCY INFORMATION C ON LUNIT. IFLAG=1 INDICATES THAT RLON AND RLAT ARE TO BE C PRINTED. C IFLAG = 1 CALL TRPRNT (N,LUNIT,RLON,RLAT,DUM,IADJ,IEND,IFLAG) C C INTERPOLATE TO AN NI BY NJ UNIFORM GRID WITH LOWER LEFT C CORNER (OLAT,OLON) AND UPPER RIGHT CORNER (OLAT+DLAT, C OLON+DLON). IFLAG=2 INDICATES THAT GRADIENT ESTIMATES C ARE TO BE COMPUTED BY UNIF AND SAVED. INTERPOLATED C VALUES ARE STORED IN WW, GRADIENT ESTIMATES IN GRAD. C CONVERT FROM DEGREES TO RADIANS. C OLAT = SC*OLAT OLON = SC*OLON DLAT = SC*DLAT DLON = SC*DLON DLAT = DLAT/FLOAT(NI-1) DLON = DLON/FLOAT(NJ-1) DO 4 I = 1,NI 4 PLAT(I) = OLAT + FLOAT(I-1)*DLAT DO 5 J = 1,NJ 5 PLON(J) = OLON + FLOAT(J-1)*DLON IFLAG = 2 CALL UNIF (N,X,Y,Z,W,IADJ,IEND,NROW,NI,NJ,PLAT, . PLON,IFLAG, GRAD,WW,IER) C C IER .GT. 0 IFF EXTRAPOLATION WAS REQUIRED C IF (IER .NE. 0) WRITE (LUNIT,101) IER C C COMPUTE INTERPOLATION ERRORS C ERMAX = 0. ERMEAN = 0. DO 6 J = 1,NJ DO 6 I = 1,NI CALL TRANS (1,PLAT(I),PLON(J), PX,PY,PZ) PW = TFUN(PX,PY,PZ) ERR = ABS(PW-WW(I,J)) ERMAX = AMAX1(ERMAX,ERR) 6 ERMEAN = ERMEAN + ERR ERMEAN = ERMEAN/FLOAT(NI*NJ) C C PRINT THE INTERPOLATION ERROR NORMS C WRITE (LUNIT,102) ERMAX, ERMEAN STOP C C FORMATS C 100 FORMAT (1H1,23HERROR IN TRMESH, IER = ,I1) 101 FORMAT (//1H ,20HERROR IN UNIF, IER =,I4) 102 FORMAT (//1H ,12HMAX ERROR = ,E10.3,15X, . 13HMEAN ERROR = ,E10.3) END