C    *******************************************************************        
C    ** FICHE F.35                                                    **        
C    ** VORONOI CONSTRUCTIONS                                         **        
C    ** TWO SEPARATE PARTS: TWO AND THREE DIMENSIONAL VERSIONS.       **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
C    *******************************************************************        
C    ** FICHE F.35  -  PART A                                         **        
C    ** THE VORONOI CONSTRUCTION IN 2D.                               **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
        PROGRAM VORON2                                                          
                                                                                
        COMMON / BLOCK1 / RX, RY                                                
                                                                                
C    *******************************************************************        
C    ** CONSTRUCTION OF THE VORONOI POLYGON IN 2D.                    **        
C    **                                                               **        
C    ** THIS PROGRAM TAKES IN A CONFIGURATION IN A SQUARE BOX WITH    **        
C    ** CONVENTIONAL PERIODIC BOUNDARY CONDITIONS AND FOR EACH ATOM   **        
C    ** OBTAINS THE SURROUNDING VORONOI POLYGON, DEFINED AS THAT      **        
C    ** REGION OF SPACE CLOSER TO THE CHOSEN ATOM THAN TO ANY OTHER.  **        
C    ** NEIGHBOURING POLYGONS DEFINE NEIGHBOURING ATOMS.              **        
C    ** THE PROGRAM IS SLOW BUT ESSENTIALLY FOOLPROOF.                **        
C    ** WE USE THE MINIMUM IMAGE CONVENTION AND SET A CUTOFF BEYOND   **        
C    ** WHICH ATOMS ARE ASSUMED NOT TO BE NEIGHBOURS: BOTH OF THESE   **        
C    ** MEASURES ARE DANGEROUS FOR SMALL AND/OR RANDOM SYSTEMS.       **        
C    ** WE DELIBERATELY DO NOT USE PREVIOUSLY-FOUND NEIGHBOURS IN     **        
C    ** CONSTRUCTING NEIGHBOUR LISTS, SO THAT AN INDEPENDENT CHECK    **        
C    ** MAY BE MADE AT THE END.                                       **        
C    ** HERE WE SIMPLY PRINT OUT THE GEOMETRICAL INFORMATION AT THE   **        
C    ** END.  THE OUTPUT IS QUITE LENGTHY.  IN PRACTICE, IT WOULD     **        
C    ** PROBABLY BE ANALYZED DIRECTLY WITHOUT PRINTING OUT.           **        
C    ** NB: BEWARE DEGENERATE CONFIGURATIONS, I.E. ONES IN WHICH MORE **        
C    ** THAN THREE VORONOI DOMAINS SHARE A VERTEX. THE SQUARE LATTICE **        
C    ** IS AN EXAMPLE.                                                **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                        NUMBER OF ATOMS              **        
C    ** REAL    RX(N),RY(N)              POSITIONS                    **        
C    ** REAL    PX(MAXCAN),PY(MAXCAN)    CANDIDATE RELATIVE POSITIONS **        
C    ** REAL    PS(MAXCAN)               SQUARED RELATIVE DISTANCES   **        
C    ** INTEGER NVER                     NUMBER OF VERTICES FOUND     **        
C    ** INTEGER NEDGE                    NUMBER OF EDGES FOUND        **        
C    ** INTEGER VERTS(MAXCAN)            VERTICES FOR EACH CANDIDATE  **        
C    **                                  = 0 IF NOT A NEIGHBOUR       **        
C    **                                  = 2 ( 1 EDGE ) IF NEIGHBOUR  **        
C    ** REAL    RXVER(MAXVER)            VERTEX RELATIVE X-COORD      **        
C    ** REAL    RYVER(MAXVER)            VERTEX RELATIVE Y-COORD      **        
C    ** INTEGER IVER(MAXVER)             ATOMIC INDICES TAGGING       **        
C    ** INTEGER JVER(MAXVER)             .. EACH VERTEX OF POLYGON    **        
C    **                                                               **        
C    ** ROUTINES REFERENCED:                                          **        
C    **                                                               **        
C    ** SUBROUTINE READCN ( CNFILE, N, BOX )                          **        
C    **    READS IN CONFIGURATION, NUMBER OF ATOMS, BOX SIZE          **        
C    ** SUBROUTINE SORT ( MAXCAN, PX, PY, PS, TAG, NCAN )             **        
C    **    SORTS NEIGHBOUR DETAILS INTO ASCENDING DISTANCE ORDER      **        
C    ** SUBROUTINE WORK ( MAXCAN, MAXVER, NCAN, NVER, NEDGE,          **        
C    **     PX, PY, PS, VERTS, RXVER, RYVER, IVER, JVER )             **        
C    **    CARRIES OUT THE VORONOI CONSTRUCTION                       **        
C    *******************************************************************        
                                                                                
        INTEGER     MAXN, MAXCAN, MAXVER                                        
        PARAMETER ( MAXN = 108, MAXCAN = 50, MAXVER = 50 )                      
                                                                                
        REAL        RX(MAXN), RY(MAXN)                                          
                                                                                
        REAL        PX(MAXCAN), PY(MAXCAN), PS(MAXCAN)                          
        INTEGER     TAG(MAXCAN), VERTS(MAXCAN)                                  
                                                                                
        REAL        RXVER(MAXVER), RYVER(MAXVER)                                
        INTEGER     IVER(MAXVER), JVER(MAXVER)                                  
        INTEGER     NABLST(MAXVER,MAXN), NNAB(MAXN), INAB, JNAB                 
                                                                                
        INTEGER     NCAN, NVER, NCOORD, NEDGE                                   
        INTEGER     I, J, CAN, VER, N                                           
        REAL        BOX, BOXINV, RCUT, RCUTSQ, COORD                            
        REAL        RXJ, RYJ, RZJ, RXIJ, RYIJ, RZIJ, RIJSQ                      
        CHARACTER   CNFILE*30                                                   
        LOGICAL     OK                                                          
                                                                                
C    *******************************************************************        
                                                                                
        WRITE(*,'(1H1,'' **** PROGRAM VORON2 ****                  '')')        
        WRITE(*,'(//1X,''VORONOI CONSTRUCTION IN 2D                '')')        
                                                                                
C    ** BASIC PARAMETERS **                                                     
                                                                                
        WRITE(*,'('' ENTER CONFIGURATION FILENAME                  '')')        
        READ (*,'(A)') CNFILE                                                   
        WRITE(*,'('' CONFIGURATION FILENAME '',A)') CNFILE                      
                                                                                
C    ** READCN MUST READ IN INITIAL CONFIGURATION  **                           
                                                                                
        CALL READCN ( CNFILE, N, BOX )                                          
                                                                                
        WRITE(*,'(1X,I5,''-ATOM CONFIGURATION'')') N                            
        WRITE(*,'('' BOX LENGTH = '',F10.5)') BOX                               
        WRITE(*,'('' ENTER NEIGHBOUR CUTOFF IN SAME UNITS '')')                 
        READ (*,*) RCUT                                                         
        WRITE(*,'('' NEIGHBOUR CUTOFF = '',F10.5)') RCUT                        
                                                                                
        RCUTSQ = RCUT ** 2                                                      
        BOXINV = 1.0 / BOX                                                      
                                                                                
C    ** ZERO ACCUMULATORS **                                                    
                                                                                
        DO 100 J = 1, N                                                         
                                                                                
           NNAB(J) = 0                                                          
                                                                                
           DO 90 INAB = 1, NVER                                                 
                                                                                
              NABLST(INAB,J) = 0                                                
                                                                                
90         CONTINUE                                                             
                                                                                
100     CONTINUE                                                                
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP STARTS                                              **        
C    *******************************************************************        
                                                                                
        DO 1000 J = 1, N                                                        
                                                                                
           IF ( MOD ( J, 2 ) .EQ. 0 ) THEN                                      
                                                                                
              WRITE(*,'(///1X,''RESULTS FOR ATOM '',I5)') J                     
                                                                                
           ELSE                                                                 
                                                                                
              WRITE(*,'(1H1,''RESULTS FOR ATOM '',I5)') J                       
                                                                                
           ENDIF                                                                
                                                                                
           RXJ = RX(J)                                                          
           RYJ = RY(J)                                                          
           CAN = 0                                                              
                                                                                
C       ** SELECT CANDIDATES **                                                 
                                                                                
           DO 500 I = 1, N                                                      
                                                                                
              IF ( I .NE. J ) THEN                                              
                                                                                
                 RXIJ = RX(I) - RXJ                                             
                 RYIJ = RY(I) - RYJ                                             
                 RXIJ = RXIJ - ANINT ( RXIJ * BOXINV ) * BOX                    
                 RYIJ = RYIJ - ANINT ( RYIJ * BOXINV ) * BOX                    
                 RIJSQ  = RXIJ ** 2 + RYIJ ** 2                                 
                                                                                
                 IF ( RIJSQ .LT. RCUTSQ ) THEN                                  
                                                                                
                    CAN = CAN + 1                                               
                                                                                
                    IF ( CAN .GT. MAXCAN ) THEN                                 
                                                                                
                       WRITE(*,'('' TOO MANY CANDIDATES '')')                   
                       STOP                                                     
                                                                                
                    ENDIF                                                       
                                                                                
                    PX(CAN)  = RXIJ                                             
                    PY(CAN)  = RYIJ                                             
                    PS(CAN)  = RIJSQ                                            
                    TAG(CAN) = I                                                
                                                                                
                 ENDIF                                                          
                                                                                
              ENDIF                                                             
                                                                                
500        CONTINUE                                                             
                                                                                
C       ** CANDIDATES HAVE BEEN SELECTED **                                     
                                                                                
           NCAN = CAN                                                           
                                                                                
C       ** SORT INTO INCREASING DISTANCE ORDER **                               
C       ** THIS SHOULD IMPROVE EFFICIENCY      **                               
                                                                                
           CALL SORT ( MAXCAN, PX, PY, PS, TAG, NCAN )                          
                                                                                
C       ** PERFORM VORONOI CONSTRUCTION **                                      
                                                                                
           CALL WORK ( MAXCAN, MAXVER, NCAN, NVER, NEDGE,                       
     :                 PX, PY, PS, VERTS,                                       
     :                 RXVER, RYVER, IVER, JVER )                               
                                                                                
C       ** WRITE OUT RESULTS **                                                 
                                                                                
           WRITE(*,'(/1X,''NUMBER OF NEIGHBOURS '',I5)') NEDGE                  
           WRITE(*,'(/1X,''NEIGHBOUR LIST '')')                                 
           WRITE(*,10001)                                                       
                                                                                
           DO 800 CAN = 1, NCAN                                                 
                                                                                
              IF ( VERTS(CAN) .NE. 0 ) THEN                                     
                                                                                
                 PS(CAN) = SQRT ( PS(CAN) )                                     
                 WRITE(*,'(1X,I5,3X,I5,3X,2F12.5,3X,F12.5)')                    
     :              TAG(CAN), VERTS(CAN), PX(CAN), PY(CAN), PS(CAN)             
                 NNAB(J) = NNAB(J) + 1                                          
                 NABLST(NNAB(J),J) = TAG(CAN)                                   
                                                                                
              ENDIF                                                             
                                                                                
800        CONTINUE                                                             
                                                                                
           WRITE(*,'(/1X,''NUMBER OF VERTICES   '',I5)') NVER                   
           WRITE(*,'(/1X,''VERTEX LIST '')')                                    
           WRITE(*,10002)                                                       
                                                                                
           DO 900 VER = 1, NVER                                                 
                                                                                
              WRITE(*,'(1X,2I5,3X,2F12.5)')                                     
     :           TAG(IVER(VER)), TAG(JVER(VER)),                                
     :           RXVER(VER), RYVER(VER)                                         
                                                                                
900        CONTINUE                                                             
                                                                                
1000    CONTINUE                                                                
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP ENDS                                                **        
C    *******************************************************************        
                                                                                
        WRITE(*,'(1H1,''FINAL SUMMARY'')')                                      
        WRITE(*,10003)                                                          
                                                                                
        NCOORD = 0                                                              
                                                                                
        DO 2000 J = 1, N                                                        
                                                                                
           NCOORD = NCOORD + NNAB(J)                                            
                                                                                
           WRITE(*,'(1X,I5,3X,I5,3X,30I3)') J, NNAB(J),                         
     :       ( NABLST(INAB,J), INAB = 1, NNAB(J) )                              
                                                                                
C       ** CHECK THAT IF I IS A NEIGHBOUR OF J **                               
C       ** THEN J IS ALSO A NEIGHBOUR OF I     **                               
                                                                                
           DO 1500 INAB = 1, NNAB(J)                                            
                                                                                
              I = NABLST(INAB,J)                                                
                                                                                
              OK = .FALSE.                                                      
              JNAB = 0                                                          
                                                                                
1200          IF ( ( .NOT. OK ) .AND. ( JNAB .LE. NNAB(I) ) ) THEN              
                                                                                
                 OK = ( J .EQ. NABLST(JNAB,I) )                                 
                 JNAB = JNAB + 1                                                
                 GOTO 1200                                                      
                                                                                
              ENDIF                                                             
                                                                                
              IF ( .NOT. OK ) THEN                                              
                                                                                
                 WRITE(*,'(1X,I3,'' IS NOT A NEIGHBOUR OF '',I3)') J, I         
                                                                                
              ENDIF                                                             
                                                                                
1500       CONTINUE                                                             
                                                                                
2000    CONTINUE                                                                
                                                                                
        COORD = REAL ( NCOORD ) / REAL ( N )                                    
                                                                                
        WRITE(*,'(/1X,'' AVERAGE COORDINATION NUMBER = '',F10.5)') COORD        
                                                                                
        STOP                                                                    
                                                                                
10001   FORMAT(/1X,'ATOM ',3X,'EDGE ',                                          
     :         /1X,'INDEX',3X,'VERTS',3X,                                       
     :         '      RELATIVE POSITION   ',3X,'  DISTANCE  ')                  
10002   FORMAT(/1X,'   INDICES         RELATIVE POSITION ')                     
10003   FORMAT(/1X,'INDEX    NABS    ... NEIGHBOUR INDICES ... ')               
                                                                                
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE READCN ( CNFILE, N, BOX )                                    
                                                                                
        COMMON / BLOCK1 / RX, RY                                                
                                                                                
C    *******************************************************************        
C    ** SUBROUTINE TO READ IN INITIAL CONFIGURATION                   **        
C    *******************************************************************        
                                                                                
        INTEGER     MAXN                                                        
        PARAMETER ( MAXN = 108 )                                                
                                                                                
        REAL        RX(MAXN), RY(MAXN), BOX                                     
        INTEGER     N                                                           
                                                                                
        CHARACTER   CNFILE*(*)                                                  
                                                                                
        INTEGER     CNUNIT, I                                                   
        PARAMETER ( CNUNIT = 10 )                                               
                                                                                
C    *******************************************************************        
                                                                                
        OPEN ( UNIT = CNUNIT, FILE = CNFILE,                                    
     :         STATUS = 'OLD', FORM = 'UNFORMATTED' )                           
                                                                                
        READ ( CNUNIT ) N, BOX                                                  
        IF ( N .GT. MAXN ) STOP ' N TOO LARGE '                                 
        READ ( CNUNIT ) ( RX(I), I = 1, N ), ( RY(I), I = 1, N )                
                                                                                
        CLOSE ( UNIT = CNUNIT )                                                 
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
                                                                                
        SUBROUTINE WORK ( MAXCAN, MAXV, NN, NV, NE, RX, RY, RS, VERTS,          
     :                    VX, VY, IV, JV )                                      
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO PERFORM VORONOI ANALYSIS                           **        
C    **                                                               **        
C    ** WE WORK INITIALLY ON DOUBLE THE CORRECT SCALE,                **        
C    ** I.E. THE EDGES OF THE POLYGON GO THROUGH THE POINTS.          **        
C    *******************************************************************        
                                                                                
        INTEGER     MAXCAN, NN, MAXV, NV, NE                                    
        INTEGER     VERTS(MAXCAN)                                               
        REAL        RX(MAXCAN), RY(MAXCAN), RS(MAXCAN)                          
        REAL        VX(MAXV), VY(MAXV)                                          
        INTEGER     IV(MAXV), JV(MAXV)                                          
                                                                                
        LOGICAL     OK                                                          
        INTEGER     I, J, L, NN1, N, V                                          
        REAL        AI, BI, CI, AJ, BJ, CJ, DET, DETINV                         
        REAL        VXIJ, VYIJ                                                  
        REAL        TOL                                                         
        PARAMETER ( TOL = 1.E-6 )                                               
                                                                                
C    *******************************************************************        
                                                                                
C    ** IF THERE ARE LESS THAN 3 POINTS GIVEN **                                
C    ** WE CANNOT CONSTRUCT A POLYGON         **                                
                                                                                
        IF ( NN .LT. 3 ) THEN                                                   
                                                                                
           WRITE(*,'('' LESS THAN 3 POINTS GIVEN TO WORK '',I5)') NN            
           STOP                                                                 
                                                                                
        ENDIF                                                                   
                                                                                
        NN1 = NN - 1                                                            
        V = 0                                                                   
                                                                                
C    ** WE AIM TO EXAMINE EACH POSSIBLE VERTEX  **                              
C    ** DEFINED BY THE INTERSECTION OF 2 EDGES  **                              
C    ** EACH EDGE IS DEFINED BY RX,RY,RS.       **                              
                                                                                
        DO 400 I = 1, NN1                                                       
                                                                                
           AI =  RX(I)                                                          
           BI =  RY(I)                                                          
           CI = -RS(I)                                                          
                                                                                
           DO 300 J = I + 1, NN                                                 
                                                                                
              AJ =  RX(J)                                                       
              BJ =  RY(J)                                                       
              CJ = -RS(J)                                                       
                                                                                
              DET = AI * BJ - AJ * BI                                           
                                                                                
              IF ( ABS ( DET ) .GT. TOL ) THEN                                  
                                                                                
C             ** THE EDGES INTERSECT **                                         
                                                                                
                 DETINV = 1.0 / DET                                             
                                                                                
                 VXIJ = ( BI * CJ - BJ * CI ) * DETINV                          
                 VYIJ = ( AJ * CI - AI * CJ ) * DETINV                          
                                                                                
C             ** NOW WE TAKE SHOTS AT THE VERTEX **                             
C             ** USING THE REMAINING EDGES ..... **                             
                                                                                
                 OK = .TRUE.                                                    
                 L  = 1                                                         
                                                                                
100              IF ( OK .AND. ( L .LE. NN ) ) THEN                             
                                                                                
                    IF ( ( L .NE. I ) .AND. ( L .NE. J ) ) THEN                 
                                                                                
                       OK = ( RX(L) * VXIJ + RY(L) * VYIJ ) .LE. RS(L)          
                                                                                
                    ENDIF                                                       
                                                                                
                    L = L + 1                                                   
                    GOTO 100                                                    
                                                                                
                 ENDIF                                                          
                                                                                
C             ** IF THE VERTEX MADE IT      **                                  
C             ** ADD IT TO THE HALL OF FAME **                                  
C             ** CONVERT TO CORRECT SCALE   **                                  
                                                                                
                 IF ( OK ) THEN                                                 
                                                                                
                    V = V + 1                                                   
                    IF ( V .GT. MAXV ) STOP 'TOO MANY VERTICES'                 
                    IV(V)  = I                                                  
                    JV(V)  = J                                                  
                    VX(V) = 0.5 * VXIJ                                          
                    VY(V) = 0.5 * VYIJ                                          
                                                                                
                 ENDIF                                                          
                                                                                
              ENDIF                                                             
                                                                                
300        CONTINUE                                                             
                                                                                
400     CONTINUE                                                                
                                                                                
C    ** THE SURVIVING VERTICES DEFINE THE VORONOI POLYGON **                    
                                                                                
        NV = V                                                                  
                                                                                
        IF ( NV .LT. 3 ) THEN                                                   
                                                                                
           WRITE(*,'('' LESS THAN 3 VERTICES FOUND IN WORK '',I5)') NV          
           STOP                                                                 
                                                                                
        ENDIF                                                                   
                                                                                
C    ** IDENTIFY NEIGHBOURING POINTS **                                         
                                                                                
        DO 500 N = 1, NN                                                        
                                                                                
           VERTS(N) = 0                                                         
                                                                                
500     CONTINUE                                                                
                                                                                
        DO 600 V = 1, NV                                                        
                                                                                
           VERTS(IV(V)) = VERTS(IV(V)) + 1                                      
           VERTS(JV(V)) = VERTS(JV(V)) + 1                                      
                                                                                
600     CONTINUE                                                                
                                                                                
C    ** POINTS WITH NONZERO VERTS ARE NEIGHBOURS **                             
C    ** IF NONZERO, VERTS SHOULD BE EQUAL TO 2   **                             
                                                                                
C    ** CHECK RESULT AND COUNT EDGES **                                         
                                                                                
        OK = .TRUE.                                                             
        NE = 0                                                                  
                                                                                
        DO 700 N = 1, NN                                                        
                                                                                
           IF ( VERTS(N) .GT. 0 ) THEN                                          
                                                                                
              NE = NE + 1                                                       
                                                                                
              IF ( VERTS(N) .NE. 2 ) THEN                                       
                                                                                
                 OK = .FALSE.                                                   
                                                                                
              ENDIF                                                             
                                                                                
           ENDIF                                                                
                                                                                
700     CONTINUE                                                                
                                                                                
        IF ( .NOT. OK ) THEN                                                    
                                                                                
           WRITE (*,'('' **** VERTEX ERROR: DEGENERACY ? **** '')')             
                                                                                
        ENDIF                                                                   
                                                                                
        IF ( NE .NE. NV ) THEN                                                  
                                                                                
           WRITE(*,'('' **** EDGE   ERROR: DEGENERACY ? ****  '')')             
                                                                                
        ENDIF                                                                   
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE SORT ( MAXCAN, RX, RY, RS, TAG, NN )                         
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO SORT NEIGHBOURS INTO INCREASING ORDER OF DISTANCE  **        
C    **                                                               **        
C    ** FOR SIMPLICITY WE USE A BUBBLE SORT - OK FOR MAXCAN SMALL.    **        
C    *******************************************************************        
                                                                                
        INTEGER MAXCAN, NN                                                      
        REAL    RX(MAXCAN), RY(MAXCAN), RS(MAXCAN)                              
        INTEGER TAG(MAXCAN)                                                     
                                                                                
        LOGICAL CHANGE                                                          
        INTEGER I, ITOP, I1, TAGI                                               
        REAL    RXI, RYI, RSI                                                   
                                                                                
C    *******************************************************************        
                                                                                
        CHANGE = .TRUE.                                                         
        ITOP = NN - 1                                                           
                                                                                
1000    IF ( CHANGE .AND. ( ITOP .GE. 1 ) ) THEN                                
                                                                                
           CHANGE = .FALSE.                                                     
                                                                                
           DO 100 I = 1, ITOP                                                   
                                                                                
              I1 = I + 1                                                        
                                                                                
              IF ( RS(I) .GT. RS(I1) ) THEN                                     
                                                                                
                 RXI = RX(I)                                                    
                 RYI = RY(I)                                                    
                 RSI = RS(I)                                                    
                 TAGI = TAG(I)                                                  
                                                                                
                 RX(I) = RX(I1)                                                 
                 RY(I) = RY(I1)                                                 
                 RS(I) = RS(I1)                                                 
                 TAG(I) = TAG(I1)                                               
                                                                                
                 RX(I1) = RXI                                                   
                 RY(I1) = RYI                                                   
                 RS(I1) = RSI                                                   
                 TAG(I1) = TAGI                                                 
                                                                                
                 CHANGE = .TRUE.                                                
                                                                                
              ENDIF                                                             
                                                                                
100        CONTINUE                                                             
                                                                                
           ITOP = ITOP - 1                                                      
           GOTO 1000                                                            
                                                                                
        ENDIF                                                                   
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
C    *******************************************************************        
C    ** FICHE F.35  - PART B                                          **        
C    ** THE VORONOI CONSTRUCTION IN 3D.                               **        
C    *******************************************************************        
                                                                                
                                                                                
                                                                                
        PROGRAM VORON3                                                          
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ                                            
                                                                                
C    *******************************************************************        
C    ** CONSTRUCTION OF VORONOI POLYHEDRON IN 3D.                     **        
C    **                                                               **        
C    ** THIS PROGRAM TAKES IN A CONFIGURATION IN A CUBIC BOX WITH     **        
C    ** CONVENTIONAL PERIODIC BOUNDARY CONDITIONS AND FOR EACH ATOM   **        
C    ** OBTAINS THE SURROUNDING VORONOI POLYHEDRON, DEFINED AS THAT   **        
C    ** REGION OF SPACE CLOSER TO THE CHOSEN ATOM THAN TO ANY OTHER.  **        
C    ** NEIGHBOURING POLYHEDRA DEFINE NEIGHBOURING ATOMS.             **        
C    ** THIS PROGRAM IS SLOW BUT ESSENTIALLY FOOLPROOF.               **        
C    ** WE USE THE MINIMUM IMAGE CONVENTION AND SET A CUTOFF BEYOND   **        
C    ** WHICH ATOMS ARE ASSUMED NOT TO BE NEIGHBOURS: BOTH OF THESE   **        
C    ** MEASURES ARE DANGEROUS FOR SMALL AND/OR RANDOM SYSTEMS.       **        
C    ** WE DELIBERATELY DO NOT USE PREVIOUSLY-FOUND NEIGHBOURS IN     **        
C    ** CONSTRUCTING NEIGHBOUR LISTS, SO THAT AN INDEPENDENT CHECK    **        
C    ** MAY BE MADE AT THE END.                                       **        
C    ** HERE WE SIMPLY PRINT OUT THE GEOMETRICAL INFORMATION AT THE   **        
C    ** END.  THE OUTPUT IS QUITE LENGTHY.  IN PRACTICE, IT WOULD     **        
C    ** PROBABLY BE ANALYZED DIRECTLY WITHOUT PRINTING IT OUT.        **        
C    ** NB: BEWARE DEGENERATE CONFIGURATIONS, I.E. ONES IN WHICH MORE **        
C    ** THAN FOUR VORONOI DOMAINS SHARE A VERTEX. THE SIMPLE CUBIC    **        
C    ** AND FACE-CENTRED CUBIC LATTICES ARE EXAMPLES.                 **        
C    **                                                               **        
C    ** PRINCIPAL VARIABLES:                                          **        
C    **                                                               **        
C    ** INTEGER N                       NUMBER OF MOLECULES           **        
C    ** REAL    RX(N),RY(N),RZ(N)       POSITIONS                     **        
C    ** REAL    PX(MAXCAN), ETC.        CANDIDATE RELATIVE POSITIONS  **        
C    ** REAL    PS(MAXCAN)              SQUARED RELATIVE DISTANCES    **        
C    ** INTEGER NVER                    NUMBER OF VERTICES FOUND      **        
C    ** INTEGER NEDGE                   NUMBER OF EDGES FOUND         **        
C    ** INTEGER NFACE                   NUMBER OF FACES FOUND         **        
C    ** INTEGER EDGES(MAXCAN)           EDGES PER FACE FOR CANDIDATES **        
C    **                                 = 0 FOR NON-NEIGHBOURS        **        
C    ** REAL    RXVER(MAXVER)           VERTEX RELATIVE X-COORD       **        
C    ** REAL    RYVER(MAXVER)           VERTEX RELATIVE Y-COORD       **        
C    ** REAL    RZVER(MAXVER)           VERTEX RELATIVE Z-COORD       **        
C    ** INTEGER IVER(MAXVER)            ATOM INDICES DEFINING         **        
C    ** INTEGER JVER(MAXVER)            .. VERTICES IN VORONOI        **        
C    ** INTEGER KVER(MAXVER)            .. POLYHEDRON.                **        
C    **                                                               **        
C    ** ROUTINES REFERENCED:                                          **        
C    **                                                               **        
C    ** SUBROUTINE READCN ( CNFILE, N, BOX )                          **        
C    **    READS CONFIGURATION, NUMBER OF ATOMS, BOX SIZE.            **        
C    ** SUBROUTINE SORT ( MAXCAN, PX, PY, PZ, PS, TAG, NCAN )         **        
C    **    SORTS NEIGHBOUR DETAILS INTO ASCENDING DISTANCE ORDER      **        
C    ** SUBROUTINE WORK ( MAXCAN, MAXVER, NCAN, NVER, NEDGE, NFACE,   **        
C    **    PX, PY, PZ, PS, EDGES,                                     **        
C    **    RXVER, RYVER, RZVER, IVER, JVER, KVER )                    **        
C    **    CARRIES OUT VORONOI CONSTRUCTION                           **        
C    *******************************************************************        
                                                                                
        INTEGER     MAXN, MAXCAN, MAXVER                                        
        PARAMETER ( MAXN = 108, MAXCAN = 50, MAXVER = 100 )                     
                                                                                
        REAL        RX(MAXN), RY(MAXN), RZ(MAXN)                                
        REAL        PX(MAXCAN), PY(MAXCAN), PZ(MAXCAN), PS(MAXCAN)              
        INTEGER     TAG(MAXCAN), EDGES(MAXCAN)                                  
        REAL        RXVER(MAXVER), RYVER(MAXVER), RZVER(MAXVER)                 
        INTEGER     IVER(MAXVER), JVER(MAXVER), KVER(MAXVER)                    
        INTEGER     NABLST(MAXVER,MAXN), NNAB(MAXN), INAB, JNAB                 
                                                                                
        INTEGER     NCAN, NVER, NEDGE, NFACE, NCOORD                            
        INTEGER     I, CAN, VER, J, N                                           
        REAL        BOX, BOXINV, RCUT, RCUTSQ, COORD                            
        REAL        RXJ, RYJ, RZJ, RXIJ, RYIJ, RZIJ, RIJSQ                      
        CHARACTER   CNFILE*30                                                   
        LOGICAL     OK                                                          
                                                                                
C    *******************************************************************        
                                                                                
        WRITE(*,'(1H1,'' **** PROGRAM VORON3 ****                  '')')        
        WRITE(*,'(//1X,''VORONOI CONSTRUCTION IN 3D                '')')        
                                                                                
C    ** BASIC PARAMETERS **                                                     
                                                                                
        WRITE(*,'('' ENTER CONFIGURATION FILENAME                  '')')        
        READ (*,'(A)') CNFILE                                                   
        WRITE(*,'('' CONFIGURATION FILENAME '',A)') CNFILE                      
                                                                                
C    ** READCN MUST READ IN INITIAL CONFIGURATION  **                           
                                                                                
        CALL READCN ( CNFILE, N, BOX )                                          
                                                                                
        WRITE(*,'(1X,I5,''-ATOM CONFIGURATION'')') N                            
        WRITE(*,'('' BOX LENGTH = '',F10.5)') BOX                               
        WRITE(*,'('' ENTER NEIGHBOUR CUTOFF IN SAME UNITS '')')                 
        READ (*,*) RCUT                                                         
        WRITE(*,'('' NEIGHBOUR CUTOFF = '',F10.5)') RCUT                        
                                                                                
        RCUTSQ = RCUT**2                                                        
        BOXINV = 1.0 / BOX                                                      
                                                                                
C    ** ZERO ACCUMULATORS **                                                    
                                                                                
        DO 100 J = 1, N                                                         
                                                                                
           NNAB(J) = 0                                                          
                                                                                
           DO 90 INAB = 1, NVER                                                 
                                                                                
              NABLST(INAB,J) = 0                                                
                                                                                
90         CONTINUE                                                             
                                                                                
100     CONTINUE                                                                
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP STARTS                                              **        
C    *******************************************************************        
                                                                                
        DO 1000 J = 1, N                                                        
                                                                                
           WRITE(*,'(1H1,'' RESULTS FOR ATOM '',I5)') J                         
                                                                                
           RXJ = RX(J)                                                          
           RYJ = RY(J)                                                          
           RZJ = RZ(J)                                                          
           CAN = 0                                                              
                                                                                
C       ** SELECT CANDIDATES **                                                 
                                                                                
           DO 500 I = 1, N                                                      
                                                                                
              IF ( I .NE. J ) THEN                                              
                                                                                
                 RXIJ = RX(I) - RXJ                                             
                 RYIJ = RY(I) - RYJ                                             
                 RZIJ = RZ(I) - RZJ                                             
                 RXIJ = RXIJ - ANINT ( RXIJ * BOXINV ) * BOX                    
                 RYIJ = RYIJ - ANINT ( RYIJ * BOXINV ) * BOX                    
                 RZIJ = RZIJ - ANINT ( RZIJ * BOXINV ) * BOX                    
                 RIJSQ  = RXIJ**2 + RYIJ**2 + RZIJ**2                           
                                                                                
                 IF ( RIJSQ .LT. RCUTSQ ) THEN                                  
                                                                                
                    CAN = CAN + 1                                               
                                                                                
                    IF ( CAN .GT. MAXCAN ) THEN                                 
                                                                                
                       WRITE(*,'('' TOO MANY CANDIDATES '')')                   
                       STOP                                                     
                                                                                
                    ENDIF                                                       
                                                                                
                    PX(CAN)  = RXIJ                                             
                    PY(CAN)  = RYIJ                                             
                    PZ(CAN)  = RZIJ                                             
                    PS(CAN)  = RIJSQ                                            
                    TAG(CAN) = I                                                
                                                                                
                 ENDIF                                                          
                                                                                
              ENDIF                                                             
                                                                                
500        CONTINUE                                                             
                                                                                
C       ** CANDIDATES HAVE BEEN SELECTED **                                     
                                                                                
           NCAN = CAN                                                           
                                                                                
C       ** SORT INTO ASCENDING ORDER OF DISTANCE **                             
C       ** THIS SHOULD IMPROVE EFFICIENCY        **                             
                                                                                
           CALL SORT ( MAXCAN, PX, PY, PZ, PS, TAG, NCAN )                      
                                                                                
C       ** PERFORM VORONOI ANALYSIS **                                          
                                                                                
           CALL WORK ( MAXCAN, MAXVER, NCAN, NVER, NEDGE, NFACE,                
     :                 PX, PY, PZ, PS, EDGES,                                   
     :                 RXVER, RYVER, RZVER, IVER, JVER, KVER )                  
                                                                                
C       ** WRITE OUT RESULTS **                                                 
                                                                                
           WRITE(*,'(/1X,''NUMBER OF NEIGHBOURS '',I5)') NFACE                  
           WRITE(*,'(/1X,''NEIGHBOUR LIST '')')                                 
           WRITE(*,10001)                                                       
                                                                                
           DO 800 CAN = 1, NCAN                                                 
                                                                                
              IF (EDGES(CAN) .NE. 0) THEN                                       
                                                                                
                 PS(CAN) = SQRT ( PS(CAN) )                                     
                 WRITE(*,'(1X,I5,3X,I5,3X,3F12.5,3X,F12.5)')                    
     :              TAG(CAN), EDGES(CAN),                                       
     :              PX(CAN), PY(CAN), PZ(CAN), PS(CAN)                          
                                                                                
                 NNAB(J) = NNAB(J) + 1                                          
                 NABLST(NNAB(J),J) = TAG(CAN)                                   
                                                                                
              ENDIF                                                             
                                                                                
800        CONTINUE                                                             
                                                                                
           WRITE(*,'(/1X,''NUMBER OF EDGES '',I5)') NEDGE                       
                                                                                
           WRITE(*,'(/1X,''NUMBER OF VERTICES '',I5)') NVER                     
           WRITE(*,'(/1X,''VERTEX LIST '')')                                    
           WRITE(*,10002)                                                       
                                                                                
           DO 900 VER = 1, NVER                                                 
                                                                                
              WRITE(*,'(1X,3I5,3X,3F12.5)')                                     
     :        TAG(IVER(VER)), TAG(JVER(VER)), TAG(KVER(VER)),                   
     :        RXVER(VER), RYVER(VER), RZVER(VER)                                
                                                                                
900        CONTINUE                                                             
                                                                                
1000    CONTINUE                                                                
                                                                                
C    *******************************************************************        
C    ** MAIN LOOP ENDS                                                **        
C    *******************************************************************        
                                                                                
        WRITE(*,'(1H1,''FINAL SUMMARY'')')                                      
        WRITE(*,10003)                                                          
                                                                                
        NCOORD = 0                                                              
                                                                                
        DO 2000 J = 1, N                                                        
                                                                                
           NCOORD = NCOORD + NNAB(J)                                            
                                                                                
           WRITE(*,'(1X,I5,3X,I5,3X,30I3)') J, NNAB(J),                         
     :        ( NABLST(INAB,J), INAB = 1, NNAB(J) )                             
                                                                                
C       ** CHECK THAT IF I IS A NEIGHBOUR OF J **                               
C       ** THEN J IS ALSO A NEIGHBOUR OF I     **                               
                                                                                
           DO 1500 INAB = 1, NNAB(J)                                            
                                                                                
              I = NABLST(INAB,J)                                                
                                                                                
              OK = .FALSE.                                                      
              JNAB = 0                                                          
                                                                                
1200          IF ( ( .NOT. OK ) .AND. ( JNAB .LE. NNAB(I) ) ) THEN              
                                                                                
                 OK = ( J .EQ. NABLST(JNAB,I) )                                 
                 JNAB = JNAB + 1                                                
                 GOTO 1200                                                      
                                                                                
              ENDIF                                                             
                                                                                
              IF ( .NOT. OK ) THEN                                              
                                                                                
                 WRITE(*,'(1X,I3,'' IS NOT A NEIGHBOUR OF '',I3)') J, I         
                                                                                
              ENDIF                                                             
                                                                                
1500       CONTINUE                                                             
                                                                                
2000    CONTINUE                                                                
                                                                                
        COORD = REAL ( NCOORD ) / REAL ( N )                                    
                                                                                
        WRITE(*,'(/1X,'' AVERAGE COORDINATION NUMBER = '',F10.5)') COORD        
                                                                                
        STOP                                                                    
                                                                                
10001   FORMAT(/1X,'ATOM ',3X,'FACE ',                                          
     :         /1X,'INDEX',3X,'EDGES',3X,                                       
     :         '            RELATIVE POSITION         ',3X,                     
     :         '  DISTANCE')                                                    
10002   FORMAT(/1X,'      INDICES           RELATIVE POSITION')                 
10003   FORMAT(/1X,'INDEX    NABS    ... NEIGHBOUR INDICES ... ')               
                                                                                
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE READCN ( CNFILE, N, BOX )                                    
                                                                                
        COMMON / BLOCK1 / RX, RY, RZ                                            
                                                                                
C    *******************************************************************        
C    ** SUBROUTINE TO READ IN INITIAL CONFIGURATION FROM UNIT 10      **        
C    *******************************************************************        
                                                                                
        INTEGER     MAXN                                                        
        PARAMETER ( MAXN = 108 )                                                
                                                                                
        REAL        RX(MAXN), RY(MAXN), RZ(MAXN), BOX                           
                                                                                
        INTEGER     CNUNIT, N, I                                                
        PARAMETER ( CNUNIT = 10 )                                               
                                                                                
        CHARACTER   CNFILE*(*)                                                  
                                                                                
                                                                                
C    *******************************************************************        
                                                                                
        OPEN ( UNIT = CNUNIT, FILE = CNFILE,                                    
     :         STATUS = 'OLD', FORM = 'UNFORMATTED' )                           
                                                                                
        READ ( CNUNIT ) N, BOX                                                  
        IF ( N .GT. MAXN ) STOP ' N TOO LARGE '                                 
        READ ( CNUNIT ) ( RX(I), I = 1, N ),                                    
     :                  ( RY(I), I = 1, N ),                                    
     :                  ( RZ(I), I = 1, N )                                     
                                                                                
        CLOSE ( UNIT = CNUNIT )                                                 
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
                                                                                
        SUBROUTINE WORK ( MAXCAN, MAXV, NN, NV, NE, NF,                         
     :                    RX, RY, RZ, RS, EDGES,                                
     :                    VX, VY, VZ, IV, JV, KV )                              
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO PERFORM VORONOI ANALYSIS                           **        
C    **                                                               **        
C    ** WE WORK INITIALLY ON DOUBLE THE CORRECT SCALE,                **        
C    ** I.E. THE FACES OF THE POLYHEDRON GO THROUGH THE POINTS.       **        
C    *******************************************************************        
                                                                                
        INTEGER     MAXCAN, NN, MAXV, NV, NE, NF                                
        INTEGER     EDGES(MAXCAN)                                               
        REAL        RX(MAXCAN), RY(MAXCAN), RZ(MAXCAN), RS(MAXCAN)              
        REAL        VX(MAXV), VY(MAXV), VZ(MAXV)                                
        INTEGER     IV(MAXV), JV(MAXV), KV(MAXV)                                
                                                                                
        LOGICAL     OK                                                          
        INTEGER     I, J, K, L, NN1, NN2, N, V                                  
        REAL        AI, BI, CI, DI, AJ, BJ, CJ, DJ, AK, BK, CK, DK              
        REAL        AB, BC, CA, DA, DB, DC, DET, DETINV                         
        REAL        VXIJK, VYIJK, VZIJK                                         
        REAL        TOL                                                         
        PARAMETER ( TOL = 1.E-6 )                                               
                                                                                
C    *******************************************************************        
                                                                                
C    ** IF THERE ARE LESS THAN 4 POINTS GIVEN **                                
C    ** WE CANNOT CONSTRUCT A POLYHEDRON      **                                
                                                                                
        IF ( NN .LT. 4 ) THEN                                                   
                                                                                
           WRITE(*,'('' LESS THAN 4 POINTS GIVEN TO WORK '',I5)') NN            
           STOP                                                                 
                                                                                
        ENDIF                                                                   
                                                                                
        NN1 = NN - 1                                                            
        NN2 = NN - 2                                                            
        V = 0                                                                   
                                                                                
C    ** WE AIM TO EXAMINE EACH POSSIBLE VERTEX  **                              
C    ** DEFINED BY THE INTERSECTION OF 3 PLANES **                              
C    ** EACH PLANE IS SPECIFIED BY RX,RY,RZ,RS  **                              
                                                                                
        DO 400 I = 1, NN2                                                       
                                                                                
           AI =  RX(I)                                                          
           BI =  RY(I)                                                          
           CI =  RZ(I)                                                          
           DI = -RS(I)                                                          
                                                                                
           DO 300 J = I + 1, NN1                                                
                                                                                
              AJ =  RX(J)                                                       
              BJ =  RY(J)                                                       
              CJ =  RZ(J)                                                       
              DJ = -RS(J)                                                       
                                                                                
              AB = AI * BJ - AJ * BI                                            
              BC = BI * CJ - BJ * CI                                            
              CA = CI * AJ - CJ * AI                                            
              DA = DI * AJ - DJ * AI                                            
              DB = DI * BJ - DJ * BI                                            
              DC = DI * CJ - DJ * CI                                            
                                                                                
              DO 200 K = J + 1, NN                                              
                                                                                
                 AK =  RX(K)                                                    
                 BK =  RY(K)                                                    
                 CK =  RZ(K)                                                    
                 DK = -RS(K)                                                    
                                                                                
                 DET = AK * BC + BK * CA + CK * AB                              
                                                                                
                 IF ( ABS ( DET ) .GT. TOL ) THEN                               
                                                                                
C                ** THE PLANES INTERSECT **                                     
                                                                                
                    DETINV = 1.0 / DET                                          
                                                                                
                    VXIJK = ( - DK * BC + BK * DC - CK * DB ) * DETINV          
                    VYIJK = ( - AK * DC - DK * CA + CK * DA ) * DETINV          
                    VZIJK = (   AK * DB - BK * DA - DK * AB ) * DETINV          
                                                                                
C                ** NOW WE TAKE SHOTS AT THE VERTEX **                          
C                ** USING THE REMAINING PLANES .... **                          
                                                                                
                    OK = .TRUE.                                                 
                    L  = 1                                                      
                                                                                
100                 IF ( OK .AND. ( L .LE. NN ) ) THEN                          
                                                                                
                       IF ( ( L .NE. I ) .AND.                                  
     :                      ( L .NE. J ) .AND.                                  
     :                      ( L .NE. K )       ) THEN                           
                                                                                
                          OK = ( ( RX(L) * VXIJK +                              
     :                             RY(L) * VYIJK +                              
     :                             RZ(L) * VZIJK  ) .LE. RS(L) )                
                                                                                
                       ENDIF                                                    
                                                                                
                       L = L + 1                                                
                       GOTO 100                                                 
                                                                                
                    ENDIF                                                       
                                                                                
C                ** IF THE VERTEX MADE IT      **                               
C                ** ADD IT TO THE HALL OF FAME **                               
C                ** CONVERT TO CORRECT SCALE   **                               
                                                                                
                    IF ( OK ) THEN                                              
                                                                                
                       V = V + 1                                                
                                                                                
                       IF ( V .GT. MAXV ) STOP 'TOO MANY VERTICES'              
                                                                                
                       IV(V)  = I                                               
                       JV(V)  = J                                               
                       KV(V)  = K                                               
                       VX(V) = 0.5 * VXIJK                                      
                       VY(V) = 0.5 * VYIJK                                      
                       VZ(V) = 0.5 * VZIJK                                      
                                                                                
                    ENDIF                                                       
                                                                                
                 ENDIF                                                          
                                                                                
200           CONTINUE                                                          
                                                                                
300        CONTINUE                                                             
                                                                                
400     CONTINUE                                                                
                                                                                
        NV = V                                                                  
                                                                                
        IF ( NV .LT. 4 ) THEN                                                   
                                                                                
           WRITE(*,'('' LESS THAN 4 VERTICES FOUND IN WORK '',I5)') NV          
           STOP                                                                 
                                                                                
        ENDIF                                                                   
                                                                                
C    ** IDENTIFY NEIGHBOURING POINTS **                                         
                                                                                
        DO 500 N = 1, NN                                                        
                                                                                
           EDGES(N) = 0                                                         
                                                                                
500     CONTINUE                                                                
                                                                                
        DO 600 V = 1, NV                                                        
                                                                                
           EDGES(IV(V)) = EDGES(IV(V)) + 1                                      
           EDGES(JV(V)) = EDGES(JV(V)) + 1                                      
           EDGES(KV(V)) = EDGES(KV(V)) + 1                                      
                                                                                
600     CONTINUE                                                                
                                                                                
C    ** POINTS WITH NONZERO EDGES ARE NEIGHBOURS **                             
                                                                                
C    ** CHECK EULER RELATION **                                                 
                                                                                
        NF = 0                                                                  
        NE = 0                                                                  
                                                                                
        DO 700 N = 1, NN                                                        
                                                                                
           IF ( EDGES(N) .GT. 0 ) NF = NF + 1                                   
           NE = NE + EDGES(N)                                                   
                                                                                
700     CONTINUE                                                                
                                                                                
        IF ( MOD ( NE, 2 ) .NE. 0 ) THEN                                        
                                                                                
           WRITE(*,'('' NONINTEGER NUMBER OF EDGES'',I5)') NE                   
           STOP                                                                 
                                                                                
        ENDIF                                                                   
                                                                                
        NE = NE / 2                                                             
                                                                                
        IF ( ( NV - NE + NF ) .NE. 2 ) THEN                                     
                                                                                
           WRITE(*,'('' **** EULER ERROR: DEGENERACY ? **** '')')               
                                                                                
        ENDIF                                                                   
                                                                                
        RETURN                                                                  
        END                                                                     
                                                                                
                                                                                
                                                                                
        SUBROUTINE SORT ( MAXCAN, RX, RY, RZ, RS, TAG, NN )                     
                                                                                
C    *******************************************************************        
C    ** ROUTINE TO SORT NEIGHBOURS INTO INCREASING ORDER OF DISTANCE  **        
C    **                                                               **        
C    ** FOR SIMPLICITY WE USE A BUBBLE SORT - OK FOR MAXCAN SMALL.    **        
C    *******************************************************************        
                                                                                
        INTEGER MAXCAN, NN                                                      
        REAL    RX(MAXCAN), RY(MAXCAN), RZ(MAXCAN), RS(MAXCAN)                  
        INTEGER TAG(MAXCAN)                                                     
                                                                                
        LOGICAL CHANGE                                                          
        INTEGER I, ITOP, I1, TAGI                                               
        REAL    RXI, RYI, RZI, RSI                                              
                                                                                
C    *******************************************************************        
                                                                                
        CHANGE = .TRUE.                                                         
        ITOP = NN - 1                                                           
                                                                                
1000    IF ( CHANGE .AND. ( ITOP .GE. 1 ) ) THEN                                
                                                                                
           CHANGE = .FALSE.                                                     
                                                                                
           DO 100 I = 1, ITOP                                                   
                                                                                
              I1 = I + 1                                                        
                                                                                
              IF ( RS(I) .GT. RS(I1) ) THEN                                     
                                                                                
                 RXI = RX(I)                                                    
                 RYI = RY(I)                                                    
                 RZI = RZ(I)                                                    
                 RSI = RS(I)                                                    
                 TAGI = TAG(I)                                                  
                                                                                
                 RX(I) = RX(I1)                                                 
                 RY(I) = RY(I1)                                                 
                 RZ(I) = RZ(I1)                                                 
                 RS(I) = RS(I1)                                                 
                 TAG(I) = TAG(I1)                                               
                                                                                
                 RX(I1) = RXI                                                   
                 RY(I1) = RYI                                                   
                 RZ(I1) = RZI                                                   
                 RS(I1) = RSI                                                   
                 TAG(I1) = TAGI                                                 
                                                                                
                 CHANGE = .TRUE.                                                
                                                                                
              ENDIF                                                             
                                                                                
100        CONTINUE                                                             
                                                                                
           ITOP = ITOP - 1                                                      
           GOTO 1000                                                            
                                                                                
        ENDIF                                                                   
                                                                                
        RETURN                                                                  
        END                                                                     
