Initial commit from Polytechnique Montreal

1 files changed, 188 insertions, 0 deletions

diff --git a/Dragon/src/PIJS3D.f b/Dragon/src/PIJS3D.f
new file mode 100644
index 0000000..9f3ded5
--- /dev/null
+++ b/Dragon/src/PIJS3D.f
@@ -0,0 +1,188 @@
+*DECK PIJS3D
+      SUBROUTINE PIJS3D(NREG,NSOUT,NSLINE,WEIGHT,RCUTOF,SIGTAL,
+     >                  SEGLEN,NRSEG,STAYIN,GOSOUT,DPR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Integration for general 3D specular tracking.
+*
+*Copyright:
+* Copyright (C) 1991 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, G. Marleau
+*
+*Parameters: input
+* NREG    total number of regions.
+* NSOUT   number of outer surface.
+* NSLINE  number of segemnts on line.
+* WEIGHT  line weight.
+* RCUTOF  MFP cut-off factor (truncate lines).
+* SIGTAL  albedo-cross section vector.
+* SEGLEN  length of track.
+* NRSEG   region crossed by track.
+*
+*Parameters: output
+* DPR     collision probabilities.
+*
+*Parameters: scratch
+* STAYIN  stay-in zone probability.
+* GOSOUT  goes-out zone probability.
+*
+*References:
+* R. Roy et al.,
+* A Cyclic Tracking Procedure for CP Calculations in 2-D Lattices
+* Conf/Advances in Math, Comp & Reactor Physics,
+* Pittsburgh, V 1, P 2.2 4-1 (1991).
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT         NONE
+*----
+* VARIABLES
+*----
+      INTEGER          NREG,NSOUT,NSLINE,NRSEG(NSLINE)
+      REAL             RCUTOF,SIGTAL(-NSOUT:NREG)
+      DOUBLE PRECISION WEIGHT,SEGLEN(NSLINE),STAYIN(NSLINE),
+     >                 GOSOUT(NSLINE)
+      DOUBLE PRECISION DPR(-NSOUT:NREG,-NSOUT:NREG)
+*----
+*  Local variables
+*----
+      INTEGER          IL,JL,NOIL,NOJL,ISODD,JSODD,IJDEL
+      DOUBLE PRECISION TTOT,XSIL,OPATH,FINV,CUTOF
+      REAL             ZERO,ONE,HALF
+      PARAMETER       (ZERO=0.0E0, ONE=1.0E0, HALF=0.5E0 )
+      REAL             SIXT,CUTEXP
+      PARAMETER       (SIXT=HALF/3.0,CUTEXP=0.02)
+      DOUBLE PRECISION EXSIL,XSIL2
+      TTOT= ONE
+*
+*1.1)    CHANGE PATHS => GOSOUT AND STAYIN PATHS, INCLUDING ALBEDOS
+*        ADD *PII* LOCAL NON-CYCLIC CONTRIBUTIONS
+      ISODD=0
+      DO 30 IL= 1, NSLINE
+        NOIL  = NRSEG(IL)
+        IF( NOIL.LT.0 )THEN
+          IF(ISODD .EQ. 1) THEN
+            ISODD=0
+*----
+*  FOR SURFACES:
+*    OLD VERSION BEFORE SURFACE DOUBLING
+*      GOSOUT= ALBEDO * SURFACE WEIGHT
+*      WHERE ALL SURFACE WEIGHTS WERE 1.0
+*    NEW VERSION WITH SURFACE DOUBLING
+*      GOSOUT= ALBEDO
+*    STAYIN = 1- ALBEDO * SURFACE WEIGHT
+*    TTOT   = PRODUCT OF GOSOUT
+*----
+            GOSOUT(IL)= SIGTAL(NOIL)
+            STAYIN(IL)= ONE - GOSOUT(IL)
+            TTOT= TTOT * GOSOUT(IL)
+          ELSE
+            ISODD=1
+*----
+*  FOR SURFACES:
+*    OLD VERSION BEFORE SURFACE DOUBLING
+*      GOSOUT= ALBEDO * SURFACE WEIGHT
+*      WHERE ALL SURFACE WEIGHTS WERE 1.0
+*    NEW VERSION WITH SURFACE DOUBLING
+*      GOSOUT= ALBEDO
+*    STAYIN = 1- ALBEDO * SURFACE WEIGHT
+*    TTOT   = PRODUCT OF GOSOUT
+*----
+            GOSOUT(IL)= SIGTAL(NOIL)
+            STAYIN(IL)= ONE
+          ENDIF
+        ELSE
+*----
+*  FOR REGIONS
+*  STAYIN = 1 -  EXP[ -CROSS SECTION * LENGTH OF NSLINE]
+*  GOSOUT = 1 -  STAYIN
+*  TTOT   = PRODUCT OF GOSOUT
+*----
+          XSIL  = SIGTAL(NOIL)
+          IF( XSIL .EQ. ZERO) THEN
+            GOSOUT(IL)= ONE
+            STAYIN(IL)= SEGLEN(IL)
+            DPR(NOIL,NOIL)= DPR(NOIL,NOIL)
+     >                         + HALF*WEIGHT*STAYIN(IL)*STAYIN(IL)
+          ELSE IF( XSIL .LT. CUTEXP) THEN
+            OPATH= SIGTAL(NOIL)*SEGLEN(IL)
+            XSIL2=OPATH*OPATH
+            EXSIL=XSIL2*(HALF-SIXT*OPATH+XSIL2/24.0)
+            STAYIN(IL)=OPATH-EXSIL
+            GOSOUT(IL)= ONE - STAYIN(IL)
+            TTOT= TTOT * GOSOUT(IL)
+            DPR(NOIL,NOIL)= DPR(NOIL,NOIL) + WEIGHT*EXSIL
+          ELSE
+            OPATH= SIGTAL(NOIL)*SEGLEN(IL)
+            EXSIL= EXP(-OPATH)
+            STAYIN(IL)= ONE - EXSIL
+            GOSOUT(IL)= EXSIL
+            TTOT= TTOT * GOSOUT(IL)
+            DPR(NOIL,NOIL)= DPR(NOIL,NOIL)
+     >                         + WEIGHT*(OPATH-STAYIN(IL))
+          ENDIF
+        ENDIF
+ 30   CONTINUE
+*
+*1.2)    COMPUTE CYCLIC FACTORS BY ANGLE
+*        USING GLOBAL TRACK ATTENUATION: BETA(TOT)*EXP(-MFP(TOT))
+      IF(TTOT .GE. ONE )THEN
+        CALL XABORT( 'PIJS3D: ALBEDOS ARE NOT COMPATIBLE')
+      ENDIF
+      FINV= WEIGHT / (ONE-TTOT)
+*
+*1.3)    ADD *PIJ* CONTRIBUTIONS FOR FORWARD SOURCES
+      ISODD=0
+      DO 50 IL= 1, NSLINE
+        NOIL  = NRSEG(IL)
+        TTOT= FINV * STAYIN(IL)
+        CUTOF= RCUTOF*TTOT
+        IF( NOIL .LT. 0) THEN
+          ISODD=MOD(ISODD+1,2)
+          JSODD=ISODD
+          DO 70 IJDEL= 1, NSLINE
+            JL= MOD(IL+IJDEL-1,NSLINE) + 1
+            NOJL=NRSEG(JL)
+            IF( NOJL .LT. 0 ) THEN
+              JSODD=MOD(JSODD+1,2)
+              IF( ISODD.EQ.1 .AND. JSODD .EQ.0) THEN
+                DPR(NOJL,NOIL)= DPR(NOJL,NOIL) + TTOT * STAYIN(JL)
+                TTOT= TTOT * GOSOUT(JL)
+                IF( TTOT.LE.CUTOF ) GO TO 55
+              ENDIF
+            ELSE IF(ISODD.EQ.1) THEN
+              DPR(NOJL,NOIL)= DPR(NOJL,NOIL) + TTOT * STAYIN(JL)
+              TTOT= TTOT * GOSOUT(JL)
+              IF( TTOT.LE.CUTOF ) GO TO 55
+            ENDIF
+ 70       CONTINUE
+        ELSE
+          JSODD=ISODD
+          DO 80 IJDEL= 1, NSLINE
+            JL= MOD(IL+IJDEL-1,NSLINE) + 1
+            NOJL=NRSEG(JL)
+            IF( NOJL .LT. 0 ) THEN
+              JSODD=MOD(JSODD+1,2)
+              IF( JSODD .EQ.0) THEN
+                DPR(NOJL,NOIL)= DPR(NOJL,NOIL) + TTOT * STAYIN(JL)
+                TTOT= TTOT * GOSOUT(JL)
+                IF( TTOT.LE.CUTOF ) GO TO 55
+              ENDIF
+            ELSE
+              DPR(NOJL,NOIL)= DPR(NOJL,NOIL) + TTOT * STAYIN(JL)
+              TTOT= TTOT * GOSOUT(JL)
+              IF( TTOT.LE.CUTOF ) GO TO 55
+            ENDIF
+ 80       CONTINUE
+        ENDIF
+ 55     CONTINUE
+ 50   CONTINUE
+      RETURN
+      END