1 files changed, 89 insertions, 0 deletions
diff --git a/Dragon/src/XELTSA.f b/Dragon/src/XELTSA.f
new file mode 100644
index 0000000..d3e32af
--- /dev/null
+++ b/Dragon/src/XELTSA.f
@@ -0,0 +1,89 @@
+*DECK XELTSA
+      SUBROUTINE XELTSA(NDIM,ITYPBC,ABSC,INDC,DENS,ANGTSA)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* To compute the integration points and periodic density
+* for cyclic tracking.
+*
+*Copyright:
+* Copyright (C) 1994 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 dimensions for the problem.
+* ITYPBC type of boundary condition (=0/2: Cartesian/hexagonal).
+* ABSC   multidimensional width of the cell.
+* INDC   index of each coordinate of the angles.
+*
+*Parameters: output
+* DENS   effective periodic density.
+* ANGTSA tracking direction and its normal.
+*
+*Reference:
+* R. Roy, G. Marleau, A. Hebert and D. Rozon,
+* A Cyclic Tracking Procedure for Collision Probability Calculations
+* in 2-D Lattice', Advances in Mathematics, Computations and 
+* Reactor Physics, Pittsburgh, PA, April 28 - May 2 (1991).
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER          NDIM,ITYPBC
+      DOUBLE PRECISION ABSC(2)
+      INTEGER          INDC(NDIM)
+      DOUBLE PRECISION DENS
+      DOUBLE PRECISION ANGTSA(NDIM,NDIM)
+*----
+*  LOCAL PARAMETERS
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='XELTSA')
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER          IDIM
+      IF(NDIM .NE. 2) CALL XABORT(NAMSBR//
+     >': Only 2-D problems permitted yet')
+      DENS=0.0D0
+      DO 10 IDIM= 1,NDIM
+         IF(ITYPBC.EQ.0) THEN
+            ! Cartesian boundary
+            ANGTSA(IDIM,1)= REAL(INDC(IDIM))/ABSC(NDIM+1-IDIM)
+         ELSE
+            ! hexagonal boundary
+            IF(IDIM.EQ.1) THEN
+              ANGTSA(IDIM,1)= REAL(INDC(IDIM))*SQRT(3.D0)
+            ELSE
+              ANGTSA(IDIM,1)= REAL(INDC(IDIM))
+            ENDIF
+         ENDIF
+         DENS= DENS + ANGTSA(IDIM,1)*ANGTSA(IDIM,1)
+   10 CONTINUE
+      DENS= SQRT(DENS)
+*----
+*  ANGTSA(*,1) is the track direction 
+*----
+      DO 20 IDIM= 1,NDIM
+         ANGTSA(IDIM,1)= ANGTSA(IDIM,1)/DENS
+   20 CONTINUE
+*----
+*  ANGTSA(*,2) is a normal to track direction
+*----
+      ANGTSA(2,2)= -ANGTSA(1,1)
+      ANGTSA(1,2)=  ANGTSA(2,1)
+*----
+*  Processing finished, return
+*----
+      RETURN
+      END