Initial commit from Polytechnique Montreal

1 files changed, 211 insertions, 0 deletions

diff --git a/Dragon/src/XELLIN.f b/Dragon/src/XELLIN.f
new file mode 100644
index 0000000..328bb67
--- /dev/null
+++ b/Dragon/src/XELLIN.f
@@ -0,0 +1,211 @@
+*DECK XELLIN
+      SUBROUTINE XELLIN(NDIM,NCP,MAXREM,REMESH,
+     >                  NSUR,NVOL,INDEL,MINDIM,MAXDIM,
+     >                  ICOORD,ICUR,INCR,TRKBEG,TRKEND,TRKDIR,
+     >                  PROJC2,TOTLEN,
+     >                  CONV,LINMAX,LENGHT,NUMERO,LINE)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Construct one tracking line which consists of two vectors.
+*
+*Copyright:
+* Copyright (C) 1987 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): R. Roy
+*
+*Parameters: input
+* NDIM    number of dimension (2 or 3).                     
+* NCP     number of cylindres of a type + 3 (< nc3max).
+* MAXREM  max number of real mesh values in REMESH.  
+* REMESH  real mesh values (rect/cyl).                 
+* NSUR    number of surfaces.                               
+* NVOL    number of zones.                                  
+* INDEL   numbering of surfaces & zones.                    
+* MINDIM  min index values for all axes (rect/cyl).    
+* MAXDIM  max index values for all axes (rect/cyl).    
+* ICOORD  principal axes direction (X/Y/Z) for meshes. 
+* ICUR    current zonal location for a track segment.  
+* INCR    increment direction for next track segment.  
+* TRKBEG  position where a track begins.               
+* TRKEND  position where a track ends.                 
+* TRKDIR  direction of a track in all axes.            
+* PROJC2  projections of TRKDIR along tracked angles.
+* TOTLEN  total lenght of the track.                   
+* CONV    segments of tracks.                          
+* LINMAX  max. number of track segments in a single track.  
+*
+*Parameters: output
+* LENGHT  relative lenght of each segment in a track.  
+* NUMERO  region identification of each track segment.  
+* LINE    lenght of the track.                          
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT          NONE
+*
+      INTEGER           NDIM, NCP, MAXREM, NSUR, NVOL, LINMAX, LINE
+      REAL              TRKBEG(NCP), TRKDIR(NCP), CONV(NCP),
+     >                  REMESH(MAXREM), TRKEND(*), PROJC2(*), TOTLEN
+      DOUBLE PRECISION  LENGHT(LINMAX)
+      INTEGER           MINDIM(NCP), MAXDIM(NCP), ICUR(NCP),
+     >                  ICOORD(NCP), INCR(NCP), INDEL(4,*),
+     >                  NUMERO(LINMAX)
+*
+      INTEGER           IORD(4), N, NP1, NP2, IBEGIN, IEND, KELVOL
+      REAL              TKBEG1, TKBEG2, TKEND1, TKEND2, R2BEG, R2END
+      DOUBLE PRECISION  CONVOK, PAT0, PAT1
+      LOGICAL           BETWEN
+      INTEGER           NEXT, NUM, I, J
+      REAL              ANORM2, CENTRE, A, B
+*
+      ANORM2(A,B)= A*A + B*B
+      CENTRE(I,J)= REMESH( MAXDIM(I-1) + J )
+      NEXT(J)= ICUR(J) + MAX( 0, INCR(J) )
+      NUM(J)= J + 1 - NSUR
+*
+*     IF THERE ARE NO CYLINDER AT ALL
+      IEND=0
+      DO 90 I   = 1, NDIM
+         N     = ICOORD(I)
+*
+*        FIND BEGINNING VOLUME #
+         ICUR(I)= MINDIM(I)
+         DO 80 J = MINDIM(I), MAXDIM(I)-1
+            IF(TRKBEG(N).GE.REMESH(J)) ICUR(I)= J
+            IF(TRKEND(N).GE.REMESH(J)) IEND= J
+   80    CONTINUE
+         IF( ICUR(I).EQ.IEND )THEN
+            CONV(I)= TOTLEN
+         ELSE
+            IF( INCR(N).EQ.0 )THEN
+               CONV(I)= TOTLEN
+            ELSE
+               CONV(I)=(REMESH(NEXT(I))-TRKBEG(N))/TRKDIR(N)
+            ENDIF
+         ENDIF
+   90 CONTINUE
+*
+      IBEGIN= MAXDIM(3) + 3
+      DO 110 I  = 4, NCP
+         N     = ICOORD(I)
+         NP1   = MOD(N  ,3) + 1
+         NP2   = MOD(N+1,3) + 1
+         TKBEG1= CENTRE(I,1) - TRKBEG(NP1)
+         TKEND1= CENTRE(I,1) - TRKEND(NP1)
+         TKBEG2= CENTRE(I,2) - TRKBEG(NP2)
+         TKEND2= CENTRE(I,2) - TRKEND(NP2)
+         R2BEG = ANORM2(TKBEG1,TKBEG2)
+         R2END = ANORM2(TKEND1,TKEND2)
+         TRKBEG(I)= (TKBEG1*TRKDIR(NP1)+TKBEG2*TRKDIR(NP2))/PROJC2(N)
+         BETWEN   =  0.0 .LT. TRKBEG(I) .AND. TRKBEG(I) .LT. TOTLEN
+         TRKDIR(I)=  R2BEG - TRKBEG(I) * TRKBEG(I) * PROJC2(N)
+         ICUR(I)  = IBEGIN - 1
+         MINDIM(I)= IBEGIN - 1
+         IEND     = IBEGIN - 1
+         DO 100 J  = IBEGIN, MAXDIM(I)
+            IF( R2BEG    .GE. REMESH(J) )ICUR(I)= J
+            IF( TRKDIR(I).GE. REMESH(J) )MINDIM(I)= J
+            IF( R2END    .GE. REMESH(J) )IEND   = J
+  100    CONTINUE
+         IBEGIN= MAXDIM(I) + 3
+         IF( ICUR(I).EQ.MINDIM(I) .AND.
+     >       IEND   .EQ.MINDIM(I)       )THEN
+            CONV(I)=TOTLEN
+         ELSE
+            IF( (BETWEN .AND. ICUR(I).NE.MINDIM(I)) .OR.
+     >                        ICUR(I).GT.IEND            )THEN
+               INCR(I)=-1
+            ELSE
+               INCR(I)=+1
+            ENDIF
+            IF( NEXT(I).GT.MAXDIM(I) )THEN
+               CONV(I)=TOTLEN
+            ELSE
+               CONV(I)= TRKBEG(I) + INCR(I) *
+     >                              SQRT((REMESH(NEXT(I))-TRKDIR(I))
+     >                            / PROJC2(ICOORD(I)))
+            ENDIF
+         ENDIF
+  110 CONTINUE
+*
+*     VOLUME TRACKED
+      LINE  = 0
+  120 LINE  = LINE + 1
+*
+*        LOOKING FOR THE MINIMUM VALUE IN "CONVOK"
+         CONVOK= TOTLEN
+         DO 130 I= 1, NCP
+            IF( CONV(I) .LT. CONVOK ) CONVOK= CONV(I)
+  130    CONTINUE
+         DO 135 I= 1, NCP
+            IORD(MIN(4,I))= ICUR(I)
+            IF(I.GT.3.AND.ICUR(I).LT.MAXDIM(I)) GOTO 136
+  135    CONTINUE
+         IORD(4)= 0
+  136    CONTINUE
+         KELVOL= NVOL
+         INDEL(1,NUM(0))= IORD(1)
+         INDEL(2,NUM(0))= IORD(2)
+         INDEL(3,NUM(0))= IORD(3)
+         INDEL(4,NUM(0))= IORD(4)
+  885    CONTINUE
+            IF( IORD(1).EQ.INDEL(1,NUM(KELVOL)).AND.
+     >          IORD(2).EQ.INDEL(2,NUM(KELVOL)).AND.
+     >          IORD(3).EQ.INDEL(3,NUM(KELVOL)).AND.
+     >          IORD(4).EQ.INDEL(4,NUM(KELVOL)) ) GO TO 895
+            KELVOL= KELVOL - 1
+         GO TO 885
+  895    IF(KELVOL.EQ.0) CALL XABORT('XELLIN: TRACKING FAILURE.')
+         NUMERO(LINE)= KELVOL
+         LENGHT(LINE)= CONVOK
+*
+*        IF "CONVOK" IS "TOTLEN" THE TRACKING IS FINISHED
+         IF( CONVOK.EQ.TOTLEN ) GO TO 160
+*
+*        UPDATE WHERE THE MINIMUM VALUE "CONVOK" IS OBTAINED
+         DO 140 I   = 1, NDIM
+            IF( CONV(I) .NE. CONVOK ) GO TO 140
+            ICUR(I)= ICUR(I) + INCR(I)
+            IF( NEXT(I).GT.MAXDIM(I) .OR.
+     >          NEXT(I).LT.MINDIM(I)      ) GO TO 160
+            N      = ICOORD(I)
+            CONV(I)= ( REMESH(NEXT(I)) - TRKBEG(N) ) / TRKDIR(N)
+  140    CONTINUE
+         DO 150 I   = 4, NCP
+            IF( CONV(I) .NE. CONVOK ) GO TO 150
+            ICUR(I)= ICUR(I) + INCR(I)
+            IF( ICUR(I) .EQ. MINDIM(I) .AND.
+     >          INCR(I) .EQ. -1              )THEN
+               INCR(I)= +1
+            ENDIF
+            IF( NEXT(I) .GT. MAXDIM(I) .OR.
+     >          NEXT(I) .LT. MINDIM(I)       )THEN
+               CONV(I)= TOTLEN
+            ELSE
+               N = ICOORD(I)
+               CONV(I)= TRKBEG(I) + INCR(I) *
+     >                  SQRT( (REMESH(NEXT(I))-TRKDIR(I))/PROJC2(N) )
+            ENDIF
+  150    CONTINUE
+*
+*     GO TO NEXT COORDINATE
+      IF( LINE .NE. LINMAX ) GO TO 120
+      CALL XABORT('XELLIN: TOO MANY TRACKS')
+  160 CONTINUE
+*
+*     TRANSFORM LOCAL COORDINATES TO PATH LENGTHS
+      PAT0=   0.0D0
+      DO 170 I= 1, LINE
+         PAT1= LENGHT(I)
+         LENGHT(I)= PAT1-PAT0
+         PAT0= PAT1
+  170 CONTINUE
+*
+      RETURN
+      END