Initial commit from Polytechnique Montreal

1 files changed, 207 insertions, 0 deletions

diff --git a/Dragon/src/XHX2D0.f b/Dragon/src/XHX2D0.f
new file mode 100644
index 0000000..844d015
--- /dev/null
+++ b/Dragon/src/XHX2D0.f
@@ -0,0 +1,207 @@
+*DECK XHX2D0
+      SUBROUTINE XHX2D0 (NGPT,ZGAUS,WGAUS,COTE,SIGT,TRONC,PII,PIS,PSS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute DP0 collision, leakage and transmission probabilities for
+* hexagonal 2D geometries.
+*
+*Copyright:
+* Copyright (C) 2002 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): M. Ouisloumen
+*
+*Parameters: input
+* NGPT    number of Gauss integration points.
+* COTE    length of one of sides of the hexagon.
+* SIGT    total cross section.
+* TRONC   voided block cutoff criterion.
+* ZGAUS   Gauss-Legendre integration points.
+* WGAUS   Gauss-Legendre integration weights.
+*
+*Parameters: output
+* PII     volume to volume reduced probability.
+* PIS     leakage probability (PIS(i) volume to side i).
+* PSS     transmission probability (PSS(i,j) side i to side j).
+*
+*Comments:
+* Faces identification for hexagon
+*                   side 4
+*                  xxxxxxxx
+*                 x        x
+*       side 5   x          x side 3
+*               x            x
+*              x              x
+*               x            x
+*       side 6   x          x side 2
+*                 x        x
+*                  xxxxxxxx
+*                   side 1
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER    NGPT
+      REAL       ZGAUS(NGPT),WGAUS(NGPT),COTE,SIGT,TRONC,PII,PIS(6),
+     +           PSS(6,6)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (MKI3=600,MKI4=600,MKI5=600)
+      PARAMETER (PI=3.141592653589793,SQRT3=1.732050807568877,
+     +           SQRT2=1.414213562373095,ALOG2=.693147180559945,
+     +           ALOG3=1.0986122886681097,ALOGX=.7676517525907618)
+*
+      REAL       TAU(3),FKI3(3),FKI4(3),FKI5(3),FKI6(3)
+      INTEGER    IROT(6,6)
+      DOUBLE PRECISION P(3),PIS10
+*----
+*  ASSUME THAT BICKLEY KI TABLES HAVE THE SAME TABULATION POINTS AND
+*  THE SAME TRUNCATION LIMIT
+*----
+      COMMON /BICKL3/BI3(0:MKI3),BI31(0:MKI3),BI32(0:MKI3),PAS3,XLIM3,L3
+      COMMON /BICKL4/BI4(0:MKI4),BI41(0:MKI4),BI42(0:MKI4),PAS4,XLIM4,L4
+      COMMON /BICKL5/BI5(0:MKI5),BI51(0:MKI5),BI52(0:MKI5),PAS5,XLIM5,L5
+*
+      SAVE IROT
+      DATA IROT /
+     +          0, 1, 2, 3, 2, 1,
+     +          1, 0, 1, 2, 3, 2,
+     +          2, 1, 0, 1, 2, 3,
+     +          3, 2, 1, 0, 1, 2,
+     +          2, 3, 2, 1, 0, 1,
+     +          1, 2, 3, 2, 1, 0/
+*
+      FUNC3(X,K)=BI3(K)+X*(BI31(K)+X*BI32(K))
+      FUNC4(X,K)=BI4(K)+X*(BI41(K)+X*BI42(K))
+      FUNC5(X,K)=BI5(K)+X*(BI51(K)+X*BI52(K))
+*----
+*  INITIALIZATION OF COLLISION PROBABILITIES
+*----
+      P(3)=0.
+      P(2)=0.
+      P(1)=0.
+*----
+*  COMPUTE CORDE =4*V*SIGMA/S=SQRT(3)*COTE*SIGMA  (AVERAGE CORDE)
+*----
+      CORDE=SQRT3*COTE*SIGT
+      IF (CORDE.LE.TRONC) GO TO 300
+*----
+*  CONSIDER EXPLICIT INTEGRATION OF F PSS
+*----
+      PI12=PI/12.
+*
+*  CONSIDER TWO CASES  1) IF ZGAUS(I)<0 => 1/COSFI>1/SINA>1/COSA
+*                      2) IF ZGAUS(I)>0 => 1/COSFI>1/COSA>1/SINA
+*
+      NGPT2=IFIX(FLOAT(NGPT)/2.)
+      DO 50 I=1,NGPT
+      FI=PI12*(1.+ZGAUS(I))
+      COSFI=COS(FI)
+      SINFI=SIN(FI)
+      COSA=SQRT3*COSFI-SINFI
+      AUX=SQRT3*SINFI
+      SINA=AUX+COSFI
+      SINB=COSFI-AUX
+*----
+*  OPTICAL LENGHTS
+*----
+      TAU(1)=CORDE/COSFI
+      IAUX1=2
+      IAUX2=3
+      IF(I.GT.NGPT2) THEN
+        IAUX1=3
+        IAUX2=2
+      ENDIF
+      TAU(IAUX1)=CORDE/SINA
+      TAU(IAUX2)=CORDE/COSA
+*
+      LB=4
+      IF(TAU(1).LT.XLIM3) THEN
+        LB=1
+      ELSEIF(TAU(2).LT.XLIM3) THEN
+        LB=2
+      ELSEIF(TAU(3).LT.XLIM3) THEN
+        LB=3
+      ENDIF
+*
+      DO 10 J=LB,3
+      K1=NINT(PAS3*TAU(J))
+      FKI3(J)=FUNC3(TAU(J),K1)
+      FKI4(J)=FUNC4(TAU(J),K1)
+      FKI5(J)=FUNC5(TAU(J),K1)
+      FKI6(J)=.8*FKI4(J)+.2*TAU(J)*(FKI3(J)-FKI5(J))
+   10 CONTINUE
+*
+      WEIGHT=0.0
+      GO TO (20,25,30,50),LB
+*----
+*  PSS CALCULATION
+*----
+   20 WEIGHT=WGAUS(I)
+      P(3)=P(3)+WEIGHT*SINB*FKI3(1)
+      P(2)=P(2)+WEIGHT*COSFI*SINA*(FKI4(IAUX1)-FKI4(1))
+      GO TO 30
+*
+   25 WEIGHT=WGAUS(I)
+      P(2)=P(2)+WEIGHT*COSFI*SINA*FKI4(IAUX1)
+   30 P(1)=P(1)+WEIGHT*SINFI*COSA*(BI4(0)-FKI4(IAUX2))
+   50 CONTINUE
+*----
+*  NORMALIZATION
+*----
+      X1=1./(3.*CORDE)
+      P(1)=X1*P(1)
+      P(2)=X1*P(2)
+      P(3)=P(3)/3.
+      PIS10=(1.-2.*(P(1)+P(2))-P(3))/(6.*CORDE)
+      PII=(1.-6.*REAL(PIS10))/SIGT
+*
+      GO TO 350
+*----
+*  USE SERIES EXPANSION FOR CORDE->0: TAYLOR SERIES OF KI FUNCTIONS
+*----
+  300 TAU0=CORDE*.5
+      TAU02=TAU0*TAU0
+      AUX=SQRT3/PI
+      AUX0=1.-.5*SQRT3
+      AUX1=AUX*ALOG3-.33333333333333333
+      AUX2=2./SQRT3-(2.+ALOGX)/3.
+      P(1)=AUX0-.5*(TAU0*AUX1-TAU02*AUX2)
+      AUX1=AUX*(2.5*ALOG3-4.*ALOG2)-.5
+      AUX2=5.*SQRT3-9.-2.*ALOG3+.5*ALOGX
+      P(2)=SQRT3-1.5+TAU0*AUX1-TAU02*AUX2/3.
+      AUX0=2.*AUX0
+      AUX1=AUX*(.5*ALOG3-ALOG2)
+      P(3)=AUX0-8.*(1./6.+AUX1)*TAU0-4.*TAU02*(AUX0-.5*ALOG3)
+      PII=COTE*SQRT3*(4.*SQRT3-8.-2.*ALOGX+10.*ALOG3)/12.
+      PIS10=(1.-SIGT*PII)/6.
+*
+  350 CONTINUE
+*----
+*  TRANSMISSION MATRIX
+*----
+      DO  59 I=1,6
+      DO  58 J=1,6
+      IB=IROT(J,I)
+      IF(IB.GT.0) THEN
+         PSS(I,J)=REAL(P(IB))
+      ELSE
+         PSS(I,J)=0.
+      ENDIF
+   58 CONTINUE
+   59 CONTINUE
+*----
+*  LEAKAGE PRABABILITIES MATRIX
+*----
+      DO 56 I=1,6
+      PIS(I)=REAL(PIS10)
+   56 CONTINUE
+*
+      RETURN
+      END