Initial commit from Polytechnique Montreal

1 files changed, 399 insertions, 0 deletions

diff --git a/Dragon/src/MCTRK.f b/Dragon/src/MCTRK.f
new file mode 100644
index 0000000..870468f
--- /dev/null
+++ b/Dragon/src/MCTRK.f
@@ -0,0 +1,399 @@
+*DECK MCTRK
+      SUBROUTINE MCTRK(IPTRK,IPRINT,NFREG,NFSUR,NDIM,NMIX,ANGBC,ITYPBC,
+     1                 MAXMSH,NUCELL,MXGSUR,MXGREG,MAXPIN,BCRT,ICODE,
+     2                 ALBEDO,IUNFLD,DGMESH,XYZL,INDEX,IDREG,DCMESH,
+     3                 ITPIN,DRAPIN,ISEED,NGRP,NL,NFM,NDEL,NED,MATCOD,
+     4                 XST,XSS,XSN2N,XSN3N,XSSNN,XSNUSI,XSCHI,XSEDI,
+     5                 MIX,ISONBR,NU,POS,ITALLY,NBSCO,NMERGE,NGCOND,
+     6                 IMERGE,INDGRP,SCORE1,SCORE2)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Simulation of a single particle from source to death.
+*
+*Copyright:
+* Copyright (C) 2008 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): B. Arsenault
+*
+*Parameters: input
+* IPTRK   pointer to the TRACKING data structure.
+* IPRINT  print level.
+* NFREG   number of regions.
+* NFSUR   number of surfaces.
+* NDIM    problem dimensions.
+* NMIX    number of mixtures in the geometry.
+* ANGBC   angular treatment for boundary conditions (=0 isotropic;
+*         =1 specular).
+* ITYPBC  type of boundary.
+* MAXMSH  maximum number of elements in MESH array.
+* NUCELL  number of cell after unfolding in $X$, $Y$ and $Z$ directions.
+* MXGSUR  maximum number of surfaces for any geometry.
+* MXGREG  maximum number of region for any geometry.
+* MAXPIN  maximum number of pins in a cell.
+* BCRT    reflection/translation array.
+* ICODE   albedo index array.
+* ALBEDO  albedo array.
+* IUNFLD  description of unfolded geometry.
+* DGMESH  meshing vector for global geometry.
+* XYZL    Cartesian boundary coordinates.
+* ISEED   the seed for the generation of random numbers.
+* NGRP    number of energy groups.
+* NL      number of Legendre orders required in the estimations
+*         (NL=1 or higher).
+* NFM     number of fissile isotopes.
+* NDEL    number of delayed precursor groups.
+* NED     number of extra edit vectors.
+* XST     total macroscopic cross sections for each mixture and energy
+*         group.
+* XSS     total scattering cross sections for each mixture and energy
+*         group.
+* XSN2N   N2N macroscopic cross sections for each mixture and energy
+*         group.
+* XSN3N   N3N macroscopic cross sections for each mixture and energy
+*         group.
+* MATCOD  region material.
+* XSSNN   in-group and out-of-group macroscopic transfert cross sections
+*         for each mixture.
+* XSNUSI  the values of Nu time the fission cross sections for each
+*         isotope per mixture and energy group.
+* XSCHI   the values of fission spectrum per isotope per mixture for
+*         each energy group.
+* XSEDI   extra edit cross sections for each mixture and energy group.
+* MIX     the mixture number where the fission occurs.
+* ISONBR  the isotopic number where the fission occurs.
+* NU      the value of the particle weight.
+* POS     location of the particle in the x, y, and z directions.
+* ITALLY  type of tally (=0 no tally; =1 score effective
+*         multiplication factor; =2 also score macrolib information).
+* NBSCO   number of macrolib related scores.
+* NMERGE  number of homogenized regions.
+* NGCOND  number of condensed energy groups.
+* IMERGE  homogenized regions indices.
+* INDGRP  condensed groups indices.
+* INDEX   undefined.
+* IDREG   undefined.
+* DCMESH  undefined.
+* ITPIN   undefined.
+* DRAPIN  undefined.
+*
+*Parameters: input/output
+* SCORE1  score for total flux and effective multiplication factor.
+* SCORE2  macrolib score matrix.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK
+      INTEGER IPRINT,NFREG,NFSUR,NDIM,NMIX,ANGBC,ITYPBC,MAXMSH,
+     1 NUCELL(3),MXGSUR,MXGREG,MAXPIN,MIX,ISONBR,
+     2 BCRT(NFSUR),ICODE(6),IUNFLD(2,NUCELL(1),NUCELL(2),NUCELL(3)),
+     3 INDEX(5,-MXGSUR:MXGREG,2),IDREG(MXGREG,2),ITPIN(3,MAXPIN),
+     4 ISEED,NGRP,NL,NFM,NDEL,NED,MATCOD(NFREG),ITALLY,NBSCO,NMERGE,
+     5 NGCOND,IMERGE(NFREG),INDGRP(NGRP)
+      REAL ALBEDO(6),XST(NMIX,NGRP),XSS(NMIX,NGRP,NL),XSN2N(NMIX,NGRP),
+     1 XSN3N(NMIX,NGRP),XSSNN(NMIX,NGRP,NGRP,NL),
+     2 XSNUSI(NMIX,NFM,NGRP,1+NDEL),XSCHI(NMIX,NFM,NGRP,1+NDEL),
+     3 XSEDI(NMIX,NGRP,NED),NU,SCORE1(3),SCORE2(NBSCO,NMERGE,NGCOND)
+      DOUBLE PRECISION DGMESH(-1:MAXMSH,4),XYZL(2,NDIM),
+     1 DCMESH(-1:MAXMSH,4,2),DRAPIN(-1:4,MAXPIN),POS(3)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER NBIND
+      PARAMETER (NBIND=7)
+      DOUBLE PRECISION XDRCST,PI,TWOPI
+      INTEGER INDX(NBIND,0:2),ODIR(3),IREG,IDIRC(2),MESHC(4,2),NSURC(2),
+     1 NREGC(2),NTPIN,IDIR,ILEV,IDS,IDSO,IPSP,IEV,II,IGROUP,IGR,
+     2 ISUR,ISUR2,IBM,IDIM,IFIRST
+      DOUBLE PRECISION CELLPO(3,2),PINCEN(3),LENGTH,MU,MUS,VIS,PHI,
+     2 VDIR(3),OMEGA(3),ALB,NUSIGF0
+      LOGICAL TRACK,ACTIVE,KILLED,VIRTUAL
+      PARAMETER (ACTIVE=.TRUE.,KILLED=.FALSE.) 
+      REAL    RAND,XSM
+      DOUBLE PRECISION PSCAT,PN2N,PN3N,XPROB,XX
+*----
+*  DATA VALUES
+*----
+      INTEGER INDOS(2,3),INDIC(-3:3)
+      DATA INDOS / 2,3,
+     1             3,1,
+     2             1,2 /
+      DATA INDIC / 5,3,1,0,2,4,6/
+*----
+*  LOAD CONSTANTS
+*----
+      PI=XDRCST('Pi',' ')
+      TWOPI=2.D0*PI
+      IFIRST=1
+*----
+*  TRACK THE NEUTRON PATH
+*----
+      DO ILEV=0,2
+         DO IDIR=1,3
+           INDX(IDIR,ILEV)=1
+         ENDDO
+         DO IDIR=4,NBIND
+           INDX(IDIR,ILEV)=0
+         ENDDO
+      ENDDO
+      ODIR(:3)=1
+
+      TRACK=ACTIVE
+      LENGTH=0.D0
+      IDSO=0
+      VIRTUAL=.FALSE.
+     
+      DO WHILE(TRACK.EQV.ACTIVE)
+        IF(IDSO.EQ.0) THEN
+          IF(LENGTH.EQ.0.D0) THEN
+            IF(MIX.EQ.-1) THEN
+              IGR=1
+            ELSE
+              CALL RANDF(ISEED,IFIRST,RAND)
+              XPROB=0.D0
+              IGR=0
+              DO WHILE((RAND.GE.XPROB).AND.(IGR.LT.NGRP))
+                IGR=IGR+1
+                XPROB=XPROB+XSCHI(MIX,ISONBR,IGR,1)
+              ENDDO
+            ENDIF
+*----
+*  SET THE TRACK DIRECTION FOR THE FIRST FLIGHT
+*----
+            XSM=0.0
+            DO II=1,NMIX
+              XSM=MAX(XSM,XST(II,IGR))
+            ENDDO
+            CALL RANDF(ISEED,IFIRST,RAND)
+            MU=2.D0*RAND-1.D0
+            MUS=DSQRT(1.D0-MU*MU)
+            CALL RANDF(ISEED,IFIRST,RAND)
+            PHI=TWOPI*RAND
+            VDIR(1)=MUS*COS(PHI)
+            VDIR(2)=MUS*SIN(PHI)
+            VDIR(3)=MU
+          ELSE
+*----
+*  IF IT IS A VIRTUAL COLLISION, KEEP THE SAME DIRECTION,
+*  IF IT IS NOT A VIRTUAL COLLISION IT CONSISTS OF AN ISOTROPIC
+*  SCATTERING REACTION.
+*----
+            IF(.NOT.VIRTUAL) THEN
+              CALL RANDF(ISEED,IFIRST,RAND)
+              PHI=TWOPI*RAND
+              CALL RANDF(ISEED,IFIRST,RAND)
+              MU=COS(PI*(2.D0*RAND-1.D0))
+              MUS=DSQRT(1.D0-MU*MU)
+              VIS=DSQRT(1.D0-VDIR(1)*VDIR(1))
+              OMEGA(1)=VDIR(1)*MU-VIS*MUS*COS(PHI)
+              OMEGA(2)=VDIR(2)*MU+MUS*(VDIR(1)*VDIR(2)*
+     1                 COS(PHI)-VDIR(3)*SIN(PHI))/VIS
+              OMEGA(3)=VDIR(3)*MU+MUS*(VDIR(1)*VDIR(3)*
+     1                 COS(PHI)+VDIR(2)*SIN(PHI))/VIS
+              DO II=1,3
+                VDIR(II)=OMEGA(II)
+              ENDDO 
+            ENDIF
+          ENDIF
+*----
+*  SAMPLE THE FREE PATH DISTANCE
+*----
+          CALL RANDF(ISEED,IFIRST,RAND)
+          LENGTH=-LOG(RAND)/XSM
+        ENDIF
+*----
+*  LOCATE THE NEUTRON IN THE GEOMETRY
+*----
+        CALL MCTPTR(IPTRK,IPRINT,NDIM,MAXMSH,ITYPBC,NUCELL,MXGSUR,
+     1     MXGREG,MAXPIN,IUNFLD,DGMESH,XYZL,NBIND,POS,LENGTH,VDIR,
+     2     ODIR,IDS,IDSO,IREG,INDX,IDIRC,MESHC,NSURC,NREGC,NTPIN,
+     3     CELLPO,PINCEN,INDEX,IDREG,DCMESH,ITPIN,DRAPIN)
+*---
+*  TALLY PROCESSING
+*---
+        IF(ITALLY.GT.0) THEN
+          CALL MCTALLY(ITALLY,NFREG,NMIX,NGRP,NL,NFM,NDEL,NED,NBSCO,
+     1    NMERGE,NGCOND,IREG,IGR,NU,MATCOD,IMERGE,INDGRP,XSM,XST,XSS,
+     2    XSN2N,XSN3N,XSSNN,XSNUSI,XSCHI,XSEDI,SCORE1,SCORE2)
+        ENDIF
+*----
+*  AN INTERACTION HAS BEEN DETECTED, DETECT IF IT IS A VIRTUAL OR A 
+*  REAL COLLISION
+*----
+        IF(IREG.GT.0) THEN
+          IBM=MATCOD(IREG)
+          CALL RANDF(ISEED,IFIRST,RAND)
+          VIRTUAL=(RAND.LE.((XSM-XST(IBM,IGR))/XSM))
+*---
+*  DETERMINE THE TYPE OF REACTION
+*---
+          IF(.NOT.VIRTUAL) THEN
+            CALL RANDF(ISEED,IFIRST,RAND)
+            PSCAT=XSS(IBM,IGR,1)/XST(IBM,IGR)
+            PN2N=XSN2N(IBM,IGR)/XST(IBM,IGR)
+            PN3N=XSN3N(IBM,IGR)/XST(IBM,IGR)
+            IF(RAND.LE.PSCAT+PN2N+PN3N) THEN
+*---
+*  ISOTROPIC SCATTERING OR NxN EVENT
+*---
+              IF(RAND.GT.PSCAT+PN2N) THEN
+                NU=3.0*NU
+              ELSE IF(RAND.GT.PSCAT) THEN
+                NU=2.0*NU
+              ENDIF
+              XX = 0.D0
+              CALL RANDF(ISEED,IFIRST,RAND)
+              IGROUP=NGRP
+              DO WHILE((RAND.GT.XX).AND.(IGROUP.GE.1))
+                XX=XX+XSSNN(IBM,IGROUP,IGR,1)/XSS(IBM,IGR,1)
+                IGROUP=IGROUP-1
+              ENDDO
+              IGR=IGROUP+1
+              TRACK=ACTIVE
+              XSM=0.0
+              DO II=1,NMIX
+                XSM=MAX(XSM,XST(II,IGR))
+              ENDDO
+            ELSE
+*---
+*  CAPTURE OR FISSION EVENT
+*---
+              TRACK = KILLED
+              MIX   = IBM
+              NUSIGF0=0.0D0
+              DO II=1,NFM
+                NUSIGF0=NUSIGF0+XSNUSI(IBM,II,IGR,1)
+              ENDDO
+              IF((NUSIGF0.GT.0.0).AND.(NFM.EQ.1)) THEN
+                NU=NU*XSNUSI(IBM,1,IGR,1)/(XST(IBM,IGR)-XSS(IBM,IGR,1)-
+     1          XSN2N(IBM,IGR)-XSN3N(IBM,IGR))
+                ISONBR=1
+              ELSE IF(NUSIGF0.GT.0.0) THEN
+                CALL RANDF(ISEED,IFIRST,RAND)
+                XPROB=0.D0
+                ISONBR=0
+                DO WHILE((RAND.GE.XPROB).AND.(ISONBR.LT.NFM))
+                  ISONBR=ISONBR+1
+                  XPROB=XPROB+XSNUSI(IBM,ISONBR,IGR,1)/NUSIGF0
+                ENDDO
+                NU=NU*XSNUSI(IBM,ISONBR,IGR,1)/(XST(IBM,IGR)-
+     1          XSS(IBM,IGR,1)-2.0*XSN2N(IBM,IGR)-3.0*XSN3N(IBM,IGR))
+              ELSE
+                NU=0.0
+                ISONBR=0
+              ENDIF
+            ENDIF
+          ENDIF
+        ELSE
+*----
+*  A BOUNDARY CONDITION HAS BEEN ENCOUNTERED
+*----
+          IF(ITYPBC.EQ.0) THEN
+*----
+*  CARTESIAN BOUNDARY
+*----
+            ISUR=-IREG
+            ISUR2=BCRT(ISUR)
+            IREG=-ISUR2
+            IF((ISUR2.EQ.ISUR).AND.(ANGBC.EQ.1)) THEN
+*             SPECULAR REFLECTIVE BOUNDARY CONDITION
+              IF(IPRINT.GT.4) WRITE(6,*) 'SPECULAR REFLECTION ON ',ISUR
+              ALB=ALBEDO(-ICODE(INDIC(IDSO)))
+              IF(ALB.EQ.0.0) THEN
+*               no leakage
+                ISONBR = 0
+                TRACK = KILLED
+              ELSE IF(ALB.EQ.1.0) THEN
+                VDIR(IDS)=-VDIR(IDS)
+              ELSE
+                CALL RANDF(ISEED,IFIRST,RAND)
+                TRACK=(RAND.LE.ALB)
+                VDIR(IDS)=-VDIR(IDS)
+              ENDIF
+            ELSE IF(ISUR2.NE.ISUR) THEN
+*             PERIODIC BOUNDARY CONDITION
+              IF(IPRINT.GT.4) WRITE(6,*) 'BC TRANSLATION FROM ',ISUR,
+     1        ' TO',ISUR2
+              IPSP=-ISUR
+              IEV=-INDIC(IDSO)
+              IF(IPRINT.GT.99) CALL MCTPSP(IPTRK,POS,IPSP,IEV)
+              IF(IDSO.GT.0) THEN
+                POS(IDS)=POS(IDS)-(XYZL(2,IDS)-XYZL(1,IDS))
+              ELSE
+                POS(IDS)=POS(IDS)+(XYZL(2,IDS)-XYZL(1,IDS))
+              ENDIF
+              IDSO=-IDSO
+            ELSE IF((ISUR2.EQ.ISUR).AND.(ANGBC.EQ.0)) THEN
+*             WHITE BOUNDARY CONDITION
+              ALB=ALBEDO(-ICODE(INDIC(IDSO)))
+              IF(ALB.EQ.0.0) THEN
+*               no leakage
+                ISONBR = 0
+                TRACK = KILLED
+              ELSE
+*               otherwise, choose randomly according to albedo
+                CALL RANDF(ISEED,IFIRST,RAND)
+                TRACK=(RAND.LE.ALB)
+                IF(TRACK.EQV.ACTIVE) THEN
+*                 FOR ISOTROPIC BC, CHOOSE RANDOMLY THE REENTERING
+*                 DIRECTION
+                  CALL RANDF(ISEED,IFIRST,RAND)
+                  MU=2.D0*RAND-1.D0
+                  MUS=DSQRT(1.D0-MU*MU)
+                  CALL RANDF(ISEED,IFIRST,RAND)
+                  PHI=PI*RAND
+                  VDIR(IDS)=-SIGN(1,IDSO)*MUS*SIN(PHI)
+                  VDIR(INDOS(1,IDS))=MUS*COS(PHI)
+                  VDIR(INDOS(2,IDS))=MU
+                  DO IDIM=1,NDIM
+                    IF(IDIM.NE.IDS) THEN
+                      CALL RANDF(ISEED,IFIRST,RAND)
+                      POS(IDIM)=XYZL(1,IDIM)+RAND*(XYZL(2,IDIM)-
+     1                XYZL(1,IDIM))
+                    ENDIF
+                  ENDDO
+                  CALL RANDF(ISEED,IFIRST,RAND)
+                  LENGTH=-LOG(RAND)/XSM
+                ELSE
+                  ISONBR = 0
+                ENDIF
+              ENDIF
+            ELSE
+              CALL XABORT('MCTRK: INVALID TYPE OF BOUNDARY CONDITION.')
+            ENDIF
+          ELSE
+*----
+*  CYLINDRICAL BOUNDARY
+*----
+            CALL XABORT('MCTRK: CYLINDRICAL BOUNDARY NOT IMPLEMENTED.')
+          ENDIF
+          IF(IPRINT.GT.99) THEN
+             IPSP=-IREG
+             IEV=-INDIC(IDSO)
+             CALL MCTPSP(IPTRK,POS,IPSP,IEV)
+          ENDIF
+        ENDIF
+*----
+*  SAVE NEUTRON PATHS IN NXT TABLE FOR PSP DISPLAY
+*----
+        IF(IPRINT.GT.99) THEN
+           IF(TRACK.EQV.KILLED) THEN
+             IPSP=-ABS(IREG)
+           ELSE
+             IPSP=ABS(IREG)
+           ENDIF
+           CALL MCTPSP(IPTRK,POS,IPSP,1)
+        ENDIF
+      ENDDO
+      RETURN
+      END