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