summaryrefslogtreecommitdiff
path: root/Dragon/src/SYBHN2.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SYBHN2.f')
-rw-r--r--Dragon/src/SYBHN2.f398
1 files changed, 398 insertions, 0 deletions
diff --git a/Dragon/src/SYBHN2.f b/Dragon/src/SYBHN2.f
new file mode 100644
index 0000000..3a7feca
--- /dev/null
+++ b/Dragon/src/SYBHN2.f
@@ -0,0 +1,398 @@
+*DECK SYBHN2
+ SUBROUTINE SYBHN2(NREG,NSURF,SIDE,Z,IZ,VOL,SIGT,TRONC,PVS,PSS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the DP-1 leakage and transmission probabilities for an
+* heterogeneous non-sectorized hexagonal cell. The tracks are computed
+* by SYBHTK.
+*
+*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): A. Hebert
+*
+*Parameters: input
+* NREG number regions in the cell.
+* NSURF number of surfaces.
+* SIDE length of one of sides of the hexagon.
+* Z real integration mesh.
+* IZ integer integration mesh.
+* VOL volumes.
+* SIGT total cross sections.
+* TRONC voided block cutoff criterion.
+*
+*Parameters: output
+* PVS leakage probability:
+* PVS(i,j) for volume i to side j with i=1,nr and j=1,18.
+* PSS transmission probability:
+* PSS(i,j) for side i to side j with i=1,18 and j=1,18.
+*
+*Reference:
+* M. Ouisloumen, Resolution par la methode des probabilites de
+* collision de l'equation integrale du transport a deux et trois
+* dimensions en geometrie hexagonale, Ph. D. thesis, Ecole
+* Polytechnique de Montreal, Montreal, October 1993.
+*
+*Comments:
+* hexagone surface identification.
+* side a,b,c
+* side 4,5,6 dir a -> isotropic
+* xxxxxxxx dir c -> tangent to surface
+* x x dir b -> normal to surface
+* side 7,8,9 x x side 1,2,3
+* x x
+* x x
+* x x
+* side 10,11,12 x x side 16,17,18
+* x x
+* xxxxxxxx
+* side 13,14,15
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER NREG,NSURF,IZ(*)
+ REAL Z(*),SIDE,SIGT(NREG),TRONC,VOL(NREG),PVS(NREG,3*NSURF),
+ 1 PSS(3*NSURF,3*NSURF)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (MKI3=600,MKI4=600,MKI5=600)
+ PARAMETER (SQ3=1.732050807568877,SQ2=1.414213562373095,
+ 1 PI=3.141592653589793,ZI30=0.785398164,ZI40=0.666666667)
+ REAL KI3,KI4,KI5
+ REAL PBB(16)
+ INTEGER IROT(18,18)
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: COSINU
+*----
+* BICKLEY TABLES
+*----
+ 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, 0, 0, 1, 2,-3, 7, 8,-9,13,14, 0, 7, 8, 9, 1, 2, 3,
+ + 0, 0, 0, 2, 4, 5, 8,10,11,14,15, 0, 8,10,-11, 2, 4,-5,
+ + 0, 0, 0, 3,-5, 6, 9,-11,12,0, 0,16,-9,11,12,-3, 5,6,
+ + 1, 2, 3, 0, 0, 0, 1, 2,-3, 7, 8,-9,13,14, 0, 7, 8, 9,
+ + 2, 4,-5, 0, 0, 0, 2, 4, 5, 8,10,11,14,15, 0, 8,10,-11,
+ + -3, 5, 6, 0, 0, 0, 3,-5, 6, 9,-11,12,0, 0,16,-9,11,12,
+ + 7, 8, 9, 1, 2, 3, 0, 0, 0, 1, 2,-3, 7, 8,-9,13,14, 0,
+ + 8,10,-11, 2, 4,-5, 0, 0, 0, 2, 4, 5, 8,10,11,14,15, 0,
+ + -9,11,12,-3, 5, 6, 0, 0, 0, 3,-5, 6, 9,-11,12,0, 0,16,
+ + 13,14, 0, 7, 8, 9, 1, 2, 3, 0, 0, 0, 1, 2,-3, 7, 8,-9,
+ + 14,15, 0, 8,10,-11, 2, 4,-5, 0, 0, 0, 2, 4, 5, 8,10,11,
+ + 0, 0,16,-9,11,12,-3, 5, 6, 0, 0, 0, 3,-5, 6, 9,-11,12,
+ + 7, 8,-9, 13,14, 0, 7, 8, 9, 1, 2, 3, 0, 0, 0, 1, 2,-3,
+ + 8,10,11,14,15, 0, 8,10,-11, 2, 4,-5, 0, 0, 0, 2, 4, 5,
+ + 9,-11,12, 0, 0,16,-9,11,12,-3, 5, 6, 0, 0, 0, 3,-5, 6,
+ + 1, 2,-3, 7, 8,-9, 13,14, 0, 7, 8, 9, 1, 2, 3, 0, 0, 0,
+ + 2, 4, 5, 8,10,11,14,15, 0, 8,10,-11, 2, 4,-5, 0, 0, 0,
+ + 3,-5, 6, 9,-11,12, 0, 0,16,-9,11,12,-3, 5, 6, 0, 0, 0/
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(COSINU(2,NSURF))
+*----
+* INTEGRATION USING THE TRACKING
+*----
+ ZERO=TRONC/(SQ3*SIDE)
+ PBB(:16)=0.0
+ PVS(:NREG,:18)=0.0
+ IZ0=2
+ IZR=4
+ DO 205 IA=1,IZ(2)
+ DO 20 I=1,NSURF
+ COSINU(1,I)=Z(IZR+1)
+ COSINU(2,I)=Z(IZR+2)
+ IZR=IZR+2
+ 20 CONTINUE
+ MNT=IZ(IZ0+1)
+ IZ0=IZ0+2
+ IZR=IZR+1
+ DO 200 IMNT=1,MNT
+ NH=IZ(IZ0+1)
+ NX=IZ(IZ0+2)
+ ISURF=IZ(IZ0+3)+1
+ JSURF=IZ(IZ0+NH+4)+1
+ DO 190 INX=1,NX
+ Z1=Z(IZR+1)/SIDE
+ IZR=IZR+1
+ DW1=Z1*COSINU(1,ISURF)
+ DW2=Z1*COSINU(1,JSURF)
+ IF((ISURF.EQ.5).AND.(JSURF.EQ.6)) THEN
+ W610=SQ3*COSINU(2,ISURF)+COSINU(1,ISURF)
+ W611=2.0*COSINU(2,JSURF)*COSINU(2,ISURF)
+ Z2=Z1*W611
+ Z3=Z1*W610
+ KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ DO 40 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+I)
+ IF(POP.LT.XLIM3) GO TO 30
+ IF(SIGTI.LE.ZERO) GO TO 50
+ PVS(III,1)=PVS(III,1)+KI3*Z1
+ PVS(III,2)=PVS(III,2)+KI4*DW1
+ GO TO 50
+ 30 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+TABKI(2,POP0)*Z(IZR+I)*Z1
+ PVS(III,2)=PVS(III,2)+KI3*Z(IZR+I)*DW1
+ ELSE
+ PVS(III,1)=PVS(III,1)+(KI3-WI3)*Z1
+ PVS(III,2)=PVS(III,2)+(KI4-WI4)*DW1
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 40 CONTINUE
+ 50 KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ K=0
+ DO 80 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+NH+1-I)
+ IF(POP.LT.XLIM3) GO TO 70
+ IF(SIGTI.LE.ZERO) GO TO 185
+ PVS(III,1)=PVS(III,1)+KI3*Z1
+ PVS(III,2)=PVS(III,2)+KI4*DW2
+ GO TO 185
+ 70 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+TABKI(2,POP0)*Z(IZR+NH+1-I)*Z1
+ PVS(III,2)=PVS(III,2)+KI3*Z(IZR+NH+1-I)*DW2
+ ELSE
+ PVS(III,1)=PVS(III,1)+(KI3-WI3)*Z1
+ PVS(III,2)=PVS(III,2)+(KI4-WI4)*DW2
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 80 CONTINUE
+ KI5=BI5(K)+POP*(BI51(K)+POP*BI52(K))
+ PBB(1)=PBB(1)+KI3*Z1
+ PBB(3)=PBB(3)+KI4*Z3
+ PBB(5)=PBB(5)+KI5*Z1
+ PBB(6)=PBB(6)+KI5*Z2
+ ELSE IF((ISURF.EQ.5).AND.(JSURF.EQ.1)) THEN
+ W610=SQ3*COSINU(1,JSURF)+COSINU(2,JSURF)
+ W511=2.0*COSINU(2,ISURF)*COSINU(2,JSURF)
+ Z2=Z1*W511
+ Z3=Z1*W610
+ KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ K=0
+ DO 100 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+I)
+ IF(POP.LT.XLIM3) GO TO 90
+ IF(SIGTI.LE.ZERO) GO TO 110
+ PVS(III,1)=PVS(III,1)+KI3*Z1
+ PVS(III,2)=PVS(III,2)+KI4*DW1
+ GO TO 110
+ 90 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+TABKI(2,POP0)*Z(IZR+I)*Z1
+ PVS(III,2)=PVS(III,2)+KI3*Z(IZR+I)*DW1
+ ELSE
+ PVS(III,1)=PVS(III,1)+(KI3-WI3)*Z1
+ PVS(III,2)=PVS(III,2)+(KI4-WI4)*DW1
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 100 CONTINUE
+ 110 KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ DO 130 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+NH+1-I)
+ IF(POP.LT.XLIM3) GO TO 120
+ IF(SIGTI.LE.ZERO) GO TO 185
+ PVS(III,1)=PVS(III,1)+KI3*Z1
+ PVS(III,2)=PVS(III,2)+KI4*DW2
+ GO TO 185
+ 120 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+TABKI(2,POP0)*Z(IZR+NH+1-I)*Z1
+ PVS(III,2)=PVS(III,2)+KI3*Z(IZR+NH+1-I)*DW2
+ ELSE
+ PVS(III,1)=PVS(III,1)+(KI3-WI3)*Z1
+ PVS(III,2)=PVS(III,2)+(KI4-WI4)*DW2
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 130 CONTINUE
+ KI5=BI5(K)+POP*(BI51(K)+POP*BI52(K))
+ PBB(7)=PBB(7)+KI3*Z1
+ PBB(9)=PBB(9)+KI4*Z3
+ PBB(11)=PBB(11)+KI5*Z1
+ PBB(12)=PBB(12)+KI5*Z2
+ ELSE IF((ISURF.EQ.5).AND.(JSURF.EQ.2)) THEN
+ Z2=Z1*COSINU(1,ISURF)
+ Z3=Z1*COSINU(1,ISURF)*COSINU(1,JSURF)
+ KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ K=0
+ DO 150 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+I)
+ IF(POP.LT.XLIM3) GO TO 140
+ IF(SIGTI.LE.ZERO) GO TO 160
+ PVS(III,1)=PVS(III,1)+KI3*Z1
+ PVS(III,2)=PVS(III,2)+KI4*Z2
+ GO TO 160
+ 140 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+TABKI(2,POP0)*Z(IZR+I)*Z1
+ PVS(III,2)=PVS(III,2)+KI3*Z(IZR+I)*Z2
+ ELSE
+ PVS(III,1)=PVS(III,1)+(KI3-WI3)*Z1
+ PVS(III,2)=PVS(III,2)+(KI4-WI4)*Z2
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 150 CONTINUE
+ 160 KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ K=0
+ DO 180 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+NH+1-I)
+ IF(POP.LT.XLIM3) GO TO 170
+ IF(SIGTI.LE.ZERO) GO TO 185
+ PVS(III,1)=PVS(III,1)+KI3*Z1
+ PVS(III,2)=PVS(III,2)+KI4*Z2
+ GO TO 185
+ 170 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+TABKI(2,POP0)*Z(IZR+NH+1-I)*Z1
+ PVS(III,2)=PVS(III,2)+KI3*Z(IZR+NH+1-I)*Z2
+ ELSE
+ PVS(III,1)=PVS(III,1)+(KI3-WI3)*Z1
+ PVS(III,2)=PVS(III,2)+(KI4-WI4)*Z2
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 180 CONTINUE
+ KI5=BI5(K)+POP*(BI51(K)+POP*BI52(K))
+ PBB(13)=PBB(13)+2.0*KI3*Z1
+ PBB(14)=PBB(14)+2.0*KI4*Z2
+ PBB(15)=PBB(15)+2.0*KI5*Z3
+ PBB(16)=PBB(16)+2.0*KI5*(Z1-Z3)
+ ENDIF
+ 185 IZR=IZR+NH
+ 190 CONTINUE
+ IZ0=IZ0+NH+4
+ 200 CONTINUE
+ 205 CONTINUE
+*----
+* PSS ORTHONORMALIZATION
+*----
+ E1=Z(1)
+ E2=Z(2)
+ E3=Z(3)
+ E4=Z(4)
+ P13=PBB(13)
+ PBB(1)=PBB(1)*E1*E1
+ PBB(7)=PBB(7)*E1*E1
+ PBB(13)=P13*E1*E1
+ PBB(3)=-0.25*PBB(3)*E1*E4*SQ3
+ PBB(9)=-0.25*PBB(9)*E1*E4
+ PBB(14)=E1*(PBB(14)*E2-E3*P13)
+ SQRT6=.25*SQ3*E4*E2
+ PBB(5)=SQRT6*PBB(5)+E3*PBB(3)/E1
+ PBB(11)=SQRT6*PBB(11)+E3*PBB(9)/E1
+ PBB(6)=-0.5*E4*E4*PBB(6)
+ PBB(12)=-0.5*E4*E4*PBB(12)
+ PBB(16)=E4*E4*PBB(16)
+ PBB(15)=E2*E2*PBB(15)-E3*(2.*PBB(14)/E1+E3*P13)
+*----
+* PIS NORMALIZATION
+*----
+ DO 210 I=1,NREG
+ COEF=0.25*SIDE/VOL(I)
+ IF(SIGT(I).LE.ZERO) THEN
+ PVS(I,2)=COEF*E1*(E2*PVS(I,2)-E3*PVS(I,1))
+ PVS(I,1)=COEF*PVS(I,1)*E1*E1
+ ELSE
+ SIGTI=COEF/SIGT(I)
+ PVS(I,2)=SIGTI*E1*(E2*PVS(I,2)-E3*PVS(I,1))
+ PVS(I,1)=SIGTI*PVS(I,1)*E1*E1
+ ENDIF
+ 210 CONTINUE
+*----
+* OTHER PROBABILITIES COMPUTATION
+*----
+ PBB(2)=(-SQ3*PBB(3)-4.*PBB(1))/SQ2
+ PBB(4)=-4.5*PBB(6)+SQ3*(-SQ2*PBB(5)+8.*PBB(3))+8.*PBB(1)
+ PBB(8)=(-3.*SQ3*PBB(9)-4.*PBB(7))/SQ2
+ PBB(10)=-4.5*PBB(12)+SQ3*(SQ2*PBB(11)+8.*PBB(9))+8.*PBB(7)
+*----
+* TRANSMISSION MATRIX
+*----
+ DO 225 I=1,18
+ DO 220 J=1,18
+ IB=IROT(J,I)
+ IF(IB.LT.0) THEN
+ PSS(I,J)=-PBB(-IB)
+ ELSEIF(IB.GT.0) THEN
+ PSS(I,J)=PBB(IB)
+ ELSE
+ PSS(I,J)=0.
+ ENDIF
+ 220 CONTINUE
+ 225 CONTINUE
+*----
+* LEAKAGE MARTIX
+*----
+ DO 235 J=1,NREG
+ K=3
+ DO 230 I=1,5
+ K=K+3
+ PVS(J,K-1)=PVS(J,2)
+ PVS(J,K-2)=PVS(J,1)
+ 230 CONTINUE
+ 235 CONTINUE
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(COSINU)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 15:16:03 +0000