diff options
Diffstat (limited to 'Dragon/src/RECT2.f')
| -rw-r--r-- | Dragon/src/RECT2.f | 371 |
1 files changed, 371 insertions, 0 deletions
diff --git a/Dragon/src/RECT2.f b/Dragon/src/RECT2.f new file mode 100644 index 0000000..4e95e43 --- /dev/null +++ b/Dragon/src/RECT2.f @@ -0,0 +1,371 @@ +*DECK RECT2 + SUBROUTINE RECT2 (NA,A,B,SIGT,TRONC,PVV,PVS,PSS,ALPA,PWA) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Integration of the collision probabilities of an homogeneous 2-D +* rectangle (DP-0 surface angular flux approximation). +* +*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): A. Hebert +* +*Parameters: input +* NA number of Gauss-Legendre base points. +* A X-thickness of the rectangle. +* B Y-thickness of the rectangle. +* SIGT cross section. +* TRONC voided block criterion. +* ALPA Gauss-Legendre base points. +* PWA Gauss-Legendre weights. +* +*Parameters: output +* PVV volume to volume reduced probability. +* PVS volume to surface probabilities: +* XINF surfaces 1, 2 and 3; XSUP surfaces 4, 5 and 6; +* YINF surfaces 7, 8 and 9; YSUP surfaces 10, 11 and 12. +* PSS surface to surface probabilities in the following order: +* PSS(i,j) is the probability from surface i to surface j. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NA + REAL A,B,SIGT,TRONC,PVV,PVS(12),PSS(12,12),ALPA(NA),PWA(NA) +*---- +* LOCAL VARIABLES +*---- + PARAMETER (MKI3=600,MKI4=600,MKI5=600) + PARAMETER (PI=3.141592654,ZI40=0.6666666667,ZI50=0.589048623, + 1 ZI60=0.533333333,COEF1=4.24264069,COEF2=2.82842712) + LOGICAL VOID + DOUBLE PRECISION Y1,Y2,AOB,BOA,R,DPREC2,DPREC5,PHI,PHJ,RM,CO,SI, + 1 Z1,Z2,POP,COSI,CO2,ZI3,ZI4,ZI5,ZI6,DEN1,DEN2,GAR(20),QUOT(2), + 2 PHIB(2),ABSC(2),PBB(28) + INTEGER ISN(12,12) + 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 ISN + DATA ISN/ 0, 0, 0, 4, 12, 0, 1, 9,-19, 1, 9, 19, + 1 0, 0, 0, 8, 16, 0, 5, 13,-23, 5, 13, 23, + 2 0, 0, 0, 0, 0, 28, 17, 21,-25,-17,-21,-25, + 3 4, 12, 0, 0, 0, 0, 1, 9, 19, 1, 9,-19, + 4 8, 16, 0, 0, 0, 0, 5, 13, 23, 5, 13,-23, + 5 0, 0, 28, 0, 0, 0,-17,-21,-25, 17, 21,-25, + 6 3, 11, 20, 3, 11,-20, 0, 0, 0, 2, 10, 0, + 7 7, 15, 24, 7, 15,-24, 0, 0, 0, 6, 14, 0, + 8 -18,-22,-27, 18, 22,-27, 0, 0, 0, 0, 0, 26, + 9 3, 11,-20, 3, 11, 20, 2, 10, 0, 0, 0, 0, + 1 7, 15,-24, 7, 15, 24, 6, 14, 0, 0, 0, 0, + 2 18, 22,-27,-18,-22,-27, 0, 0, 26, 0, 0, 0/ +* + AOB=A/B + BOA=B/A + X=2.0*A*B/(A+B) + VOID=(X*SIGT.LE.TRONC) + SIG=1.0E10 + IF(VOID) GO TO 105 + SIG=1.0/SIGT + IF(A.EQ.B) GO TO 60 +*---- +* RECTANGULAR CELL +*---- + QUOT(1)=BOA + QUOT(2)=AOB + PHIB(1)=ATAN(AOB) + PHIB(2)=ATAN(BOA) + ABSC(1)=B + ABSC(2)=A + DO 10 I=1,20 + GAR(I)=0.0 +10 CONTINUE + DO 55 IL=1,NA + X=0.5*(ALPA(IL)+1.0) + WA=PWA(IL) + L=0 + DO 50 IZ=1,2 + PHI=PHIB(IZ)*X + CO=COS(PHI) + SI=SIN(PHI) + COSI=CO*SI + CO2=CO*CO + Z1=PHIB(IZ)*WA*COSI/(PI*A) + Z2=PHIB(IZ)*WA*(CO-QUOT(IZ)*SI)*2.0/PI + POP=SIGT*ABSC(IZ)/CO + ZI3=0.0 + ZI4=0.0 + ZI5=0.0 + IF(POP.GE.XLIM3) GO TO 40 + K=NINT(POP*PAS3) + ZI3=BI3(K)+POP*(BI31(K)+POP*BI32(K)) + ZI4=BI4(K)+POP*(BI41(K)+POP*BI42(K)) + ZI5=BI5(K)+POP*(BI51(K)+POP*BI52(K)) +40 ZI6=0.8*ZI4+0.2*POP*(ZI3-ZI5) + GAR(L+1)=GAR(L+1)+(ZI40-ZI4)*Z1 + GAR(L+2)=GAR(L+2)+(ZI50-ZI5)*Z1*CO + GAR(L+3)=GAR(L+3)+(ZI50-ZI5)*Z1*SI + GAR(L+4)=GAR(L+4)+(ZI60-ZI6)*Z1*COSI + GAR(L+5)=GAR(L+5)+(ZI60-ZI6)*Z1*CO2 + GAR(L+6)=GAR(L+6)+(ZI60-ZI6)*Z1 + GAR(L+7)=GAR(L+7)+ZI3*Z2 + GAR(L+8)=GAR(L+8)+ZI4*Z2*CO + GAR(L+9)=GAR(L+9)+ZI5*Z2*CO2 + GAR(L+10)=GAR(L+10)+ZI5*Z2 + L=L+10 +50 CONTINUE +55 CONTINUE + PBB(1)=(GAR(1)+GAR(11))*SIG + PBB(2)=GAR(7) + PBB(4)=GAR(17) + PBB(5)=(GAR(3)+GAR(12))*SIG + PBB(6)=GAR(8) + PBB(8)=GAR(18) + PBB(9)=(GAR(2)+GAR(13))*SIG + PBB(10)=GAR(8) + PBB(12)=GAR(18) + PBB(13)=(GAR(4)+GAR(14))*SIG + PBB(14)=GAR(9) + PBB(16)=GAR(19) + PBB(17)=PBB(9) + PBB(19)=PBB(5) + PBB(21)=(GAR(5)+GAR(16)-GAR(15))*SIG + PBB(23)=(GAR(6)-GAR(5)+GAR(15))*SIG + PBB(25)=PBB(13) + PBB(26)=GAR(10)-GAR(9) + PBB(28)=GAR(20)-GAR(19) + GO TO 110 +*---- +* SQUARE CELL +*---- +60 PHI1=0.25*PI + DO 65 I=1,8 + GAR(I)=0.0 +65 CONTINUE + DO 100 IL=1,NA + X=0.5*(ALPA(IL)+1.0) + WA=PWA(IL) + PHI=PHI1*X + CO=COS(PHI) + SI=SIN(PHI) + COSI=CO*SI + Z1=0.25*WA*COSI/A + Z2=0.5*WA*(CO-SI) + POP=SIGT*A/CO + ZI3=0.0 + ZI4=0.0 + ZI5=0.0 + IF(POP.GE.XLIM3) GO TO 90 + K=NINT(POP*PAS3) + ZI3=BI3(K)+POP*(BI31(K)+POP*BI32(K)) + ZI4=BI4(K)+POP*(BI41(K)+POP*BI42(K)) + ZI5=BI5(K)+POP*(BI51(K)+POP*BI52(K)) +90 ZI6=0.8*ZI4+0.2*POP*(ZI3-ZI5) + GAR(1)=GAR(1)+(ZI40-ZI4)*Z1 + GAR(2)=GAR(2)+(ZI50-ZI5)*Z1*(CO+SI) + GAR(3)=GAR(3)+(ZI60-ZI6)*Z1*COSI + GAR(4)=GAR(4)+(ZI60-ZI6)*Z1 + GAR(5)=GAR(5)+ZI3*Z2 + GAR(6)=GAR(6)+ZI4*Z2*CO + GAR(7)=GAR(7)+ZI5*Z2*CO*CO + GAR(8)=GAR(8)+ZI5*Z2 +100 CONTINUE + PBB(1)=2.0*GAR(1)*SIG + PBB(2)=GAR(5) + PBB(4)=PBB(2) + PBB(5)=GAR(2)*SIG + PBB(6)=GAR(6) + PBB(8)=PBB(6) + PBB(9)=PBB(5) + PBB(10)=PBB(6) + PBB(12)=PBB(6) + PBB(13)=2.0*GAR(3)*SIG + PBB(14)=GAR(7) + PBB(16)=PBB(14) + PBB(17)=GAR(2)*SIG + PBB(19)=PBB(17) + PBB(21)=GAR(4)*SIG + PBB(23)=PBB(21) + PBB(25)=2.0*GAR(3)*SIG + PBB(26)=GAR(8)-GAR(7) + PBB(28)=PBB(26) + GO TO 110 +*---- +* VOIDED CELL +*---- +105 Z1=1.0/(8.0*A*B) + IF(A.GT.1.0E2*B) THEN + PHI=0.5D0*PI + R=A*(1.0D0+BOA*BOA*(0.5D0-0.125D0*BOA*BOA)) + RM=BOA*(0.5D0-0.125D0*BOA*BOA) + Y1=LOG(DBLE((R+A)/B)) + PBB(1)=(0.5D0*(1.0D0-RM)+SIGT*B*LOG(B/R)/PI)*BOA + PBB(2)=(R-B)/A-4.0D0*SIGT*(B*PHI+B*BOA*LOG(B/R))/PI + 1 +SIGT*SIGT*B*B*(Y1-(R-B)/A) + PBB(4)=RM + PBB(5)=((1.0-RM)/3.0-0.25*SIGT*B*Y1)*BOA + PBB(6)=4.0*PHI/(3.0*PI)-SIGT*BOA*(R-B) + 1 +2.0*SIGT*SIGT*(B*B*PHI+B*B*BOA*LOG(B/R))/PI + PBB(8)=RM*2.0/3.0 + PBB(9)=(-0.25*SIGT*B*(1.0-BOA*2.0/3.0))*BOA + PBB(13)=(0.125*A/R-SIGT*B/6.0)*BOA + PBB(14)=(R*(R-B)+A*A)/(4.0*A*R)-4.0*SIGT*B*PHI/(3.0*PI) + 1 +0.5*SIGT*SIGT*B*BOA*(R-B) + PBB(16)=RM*4.0/9.0 + PBB(21)=0.125*(1.0-RM*RM*B/R)*BOA + PBB(23)=0.125*((B+B-B*B/R)/B-RM)*BOA + PBB(26)=0.25*PBB(2)*(R-B)/R + PBB(28)=0.25*RM*RM*B/R + PVS(1)=REAL(-0.5D0*BOA*LOG(B/R)/PI) + PVS(7)=REAL((PHI+0.5D0*BOA*LOG(B/R))/PI-0.25D0*SIGT*(B*Y1- + 1 BOA*(R-B))) + DPREC2=Z1*B*(2.0*(R-B)+B*(1.0-BOA*2.0/3.0)) + 1 -0.5*SIGT*B*(PHI+BOA*LOG(B/R))/PI + DPREC5=Z1*B*B*Y1 + PVV=REAL(0.5D0*B*Y1-0.5D0*BOA*(R-B)) + ELSE IF(B.GT.1.0E2*A) THEN + PHJ=0.5D0*PI + R=B*(1.0D0+AOB*AOB*(0.5D0-0.125D0*AOB*AOB)) + RM=AOB*(0.5D0-0.125D0*AOB*AOB) + Y2=LOG(DBLE((R+B)/A)) + PBB(1)=0.5D0*(1.0D0-RM)+SIGT*A*LOG(A/R)/PI + PBB(2)=RM + PBB(4)=(R-A)/B-4.0D0*SIGT*(A*PHJ+A*AOB*LOG(A/R))/PI + 1 +SIGT*SIGT*A*A*(Y2-(R-A)/B) + PBB(5)=-0.25*SIGT*A*(1.0-AOB*2.0/3.0) + PBB(6)=RM*2.0/3.0 + PBB(8)=4.0*PHJ/(3.0*PI)-SIGT*AOB*(R-A) + 1 +2.0*SIGT*SIGT*(A*A*PHJ+A*A*AOB*LOG(A/R))/PI + PBB(9)=(1.0-RM)/3.0-0.25*SIGT*A*Y2 + PBB(13)=0.125*B/R-SIGT*A/6.0 + PBB(14)=RM*4.0/9.0 + PBB(16)=(R*(R-A)+B*B)/(4.0*B*R)-4.0*SIGT*A*PHJ/(3.0*PI) + 1 +0.5*SIGT*SIGT*A*AOB*(R-A) + PBB(21)=0.125*((A+A-A*A/R)/A-RM) + PBB(23)=0.125*(1.0-RM*RM*A/R) + PBB(26)=0.25*RM*RM*A/R + PBB(28)=0.25*PBB(4)*(R-A)/R + PVS(1)=REAL((PHJ+0.5D0*AOB*LOG(A/R))/PI-0.25D0*SIGT*(A*Y2- + 1 AOB*(R-A))) + PVS(7)=REAL(-0.5D0*AOB*LOG(A/R)/PI) + DPREC2=Z1*A*A*Y2 + DPREC5=Z1*A*(2.0*(R-A)+A*(1.0-AOB*2.0/3.0)) + 1 -0.5*SIGT*A*(PHJ+AOB*LOG(A/R))/PI + PVV=REAL(0.5D0*A*Y2-0.5D0*AOB*(R-A)) + ELSE + PHI=ATAN(AOB) + PHJ=ATAN(BOA) + R=SQRT(A*A+B*B) + Y1=LOG(DBLE((R+A)/B)) + Y2=LOG(DBLE((R+B)/A)) + PBB(1)=(A+B-R)/(A+A)+SIGT*(A*LOG(A/R)+B*BOA*LOG(B/R))/PI + 1 +SIGT*SIGT*(R*R*R-A*A*A-B*B*B)/(6.0*A) + PBB(2)=(R-B)/A-4.0*SIGT*(B*PHI+B*BOA*LOG(B/R))/PI + 1 +SIGT*SIGT*B*B*(Y1-(R-B)/A) + PBB(4)=(R-A)/B-4.0*SIGT*A*(PHJ+AOB*LOG(A/R))/PI + 1 +SIGT*SIGT*A*A*(Y2-(R-A)/B) + PBB(5)=(2.0*PHI*B)/(3.0*PI*A)+0.25*SIGT*(-B*BOA*Y1+R-A) + 1 +SIGT*SIGT*(B*B-B*B*BOA*PHI-A*A*LOG(A/R))/(3.0*PI) + PBB(6)=4.0*PHI/(3.0*PI)-SIGT*BOA*(R-B) + 1 +2.0*SIGT*SIGT*(B*B*PHI+B*B*BOA*LOG(B/R))/PI + PBB(8)=4.0*PHJ/(3.0*PI)-SIGT*AOB*(R-A) + 1 +2.0*SIGT*SIGT*A*A*(PHJ+AOB*LOG(A/R))/PI + PBB(9)=2.0*PHJ/(3.0*PI)+0.25*SIGT*(-A*Y2+BOA*(R-B)) + 1 +SIGT*SIGT*(-B*B*BOA*LOG(B/R)+A*B-A*A*PHJ)/(3.0*PI) + PBB(13)=B/(8.0*R)-SIGT*(B*BOA*PHI+A*PHJ-B)/(3.0*PI) + 1 +SIGT*SIGT*(B*B*BOA*Y1+A*A*Y2-B*R)/12.0 + PBB(14)=(R*(R-B)+A*A)/(4.0*A*R)-4.0*SIGT*B*PHI/(3.0*PI) + 1 +0.5*SIGT*SIGT*B*BOA*(R-B) + PBB(16)=(R*(R-A)+B*B)/(4.0*B*R)-4.0*SIGT*A*PHJ/(3.0*PI) + 1 +0.5*SIGT*SIGT*A*AOB*(R-A) + PBB(21)=(A+A+B-A*A/R-R)/(8.0*A)+2.0*SIGT*A*LOG(A/R)/(3.0*PI) + 1 +SIGT*SIGT*(2.0*A*(R-A)+B*BOA*(R-B))/12.0 + PBB(23)=(A+B+B-B*B/R-R)/(8.0*A)+2.0*SIGT*B*BOA*LOG(B/R)/ + 1 (3.0*PI)+SIGT*SIGT*(2.0*B*BOA*(R-B)-A*(R-A))/12.0 + PBB(26)=0.25*PBB(2)*(R-B)/R + PBB(28)=0.25*PBB(4)*(R-A)/R + PVS(1)=REAL((0.5*AOB*LOG(A/R)-0.5*BOA*LOG(B/R)+PHJ)/PI + 1 -0.25*SIGT*((R*R*R-A*A*A-B*B*B)/(3.0*A*B)-(R-A)*AOB+A*Y2)) + PVS(7)=REAL((0.5*BOA*LOG(B/R)-0.5*AOB*LOG(A/R)+PHI)/PI + 1 -0.25*SIGT*((R*R*R-A*A*A-B*B*B)/(3.0*A*B)-(R-B)*BOA+B*Y1)) + DPREC2=Z1*(B*(R-B)+A*A*Y2)-SIGT*(A+2.0*B*BOA*LOG(B/R) + 1 +3.0*B*PHI-A*AOB*PHJ)/(6.0*PI) + DPREC5=Z1*(A*(R-A)+B*B*Y1)-SIGT*(B+2.0*A*AOB*LOG(A/R) + 1 +3.0*A*PHJ-B*BOA*PHI)/(6.0*PI) + PVV=REAL((A*A*A+B*B*B-R*R*R)/(6.0*A*B)+0.5*B*Y1+0.5*A*Y2) + ENDIF + PBB(10)=PBB(6) + PBB(12)=PBB(8) + PBB(17)=PBB(9) + PBB(19)=PBB(5) + PBB(25)=PBB(13) + PVS(2)=REAL(COEF1*DPREC5-COEF2*PVS(1)) + PVS(3)=0.0 + PVS(8)=REAL(COEF1*DPREC2-COEF2*PVS(7)) + PVS(9)=0.0 +* +110 PBB(3)=PBB(1)*AOB + PBB(7)=PBB(9)*AOB + PBB(11)=PBB(5)*AOB + PBB(15)=PBB(13)*AOB + PBB(18)=PBB(19)*AOB + PBB(20)=PBB(17)*AOB + PBB(22)=PBB(23)*AOB + PBB(24)=PBB(21)*AOB + PBB(27)=PBB(25)*AOB +*---- +* ORTHONORMALIZATION +*---- + DO 120 I=1,4 + DEN1=PBB(4+I) + DEN2=PBB(8+I) + PBB(4+I)=COEF1*DEN1-COEF2*PBB(I) + PBB(8+I)=COEF1*DEN2-COEF2*PBB(I) + PBB(12+I)=18.0*PBB(12+I)-12.0*(DEN1+DEN2)+8.0*PBB(I) + PBB(24+I)=4.0*PBB(24+I) +120 CONTINUE + DO 130 I=1,2 + DEN1=PBB(16+I) + DEN2=PBB(18+I) + PBB(16+I)=2.0*DEN1 + PBB(18+I)=2.0*DEN2 + PBB(20+I)=2.0*(COEF1*PBB(20+I)-COEF2*DEN1) + PBB(22+I)=2.0*(COEF1*PBB(22+I)-COEF2*DEN2) +130 CONTINUE +* + IF(.NOT.VOID) THEN + PVS(7)=REAL(0.25*SIG*(1.0-2.0*PBB(1)-PBB(2))/B) + PVS(8)=REAL(-0.25*SIG*(2.0*PBB(9)+PBB(10))/B) + PVS(9)=0.0 + PVS(1)=REAL(0.25*SIG*(1.0-2.0*PBB(3)-PBB(4))/A) + PVS(2)=REAL(-0.25*SIG*(2.0*PBB(11)+PBB(12))/A) + PVS(3)=0.0 + PVV=REAL(SIG*(1.0-2.0*(PVS(7)+PVS(1)))) + ENDIF + PVS(4)=PVS(1) + PVS(5)=PVS(2) + PVS(6)=PVS(3) + PVS(10)=PVS(7) + PVS(11)=PVS(8) + PVS(12)=PVS(9) + DO 150 JC=1,12 + DO 140 IC=1,12 + PSS(IC,JC)=0.0 + IB=ISN(IC,JC) + IF(IB.LT.0) THEN + PSS(IC,JC)=REAL(-PBB(-IB)) + ELSE IF(IB.GT.0) THEN + PSS(IC,JC)=REAL(PBB(IB)) + ENDIF + 140 CONTINUE + 150 CONTINUE + RETURN + END |