From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/XHX2D1.f | 322 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 322 insertions(+) create mode 100644 Dragon/src/XHX2D1.f (limited to 'Dragon/src/XHX2D1.f') diff --git a/Dragon/src/XHX2D1.f b/Dragon/src/XHX2D1.f new file mode 100644 index 0000000..065e1c1 --- /dev/null +++ b/Dragon/src/XHX2D1.f @@ -0,0 +1,322 @@ +*DECK XHX2D1 + SUBROUTINE XHX2D1 (NGPT,ZGAUS,WGAUS,COTE,SIGT,TRONC,PII,PIS,PSS, + + P) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute DP1 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). +* +*Parameters: scratch +* P undefined. +* +*Comments: +* Faces identification for hexagon +* 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 +*---- + REAL FUNC3,FUNC4,FUNC5,X + INTEGER NGPT + REAL ZGAUS(NGPT),WGAUS(NGPT),COTE,SIGT,TRONC,PII,PIS(18), + + PSS(18,18) + DOUBLE PRECISION P(16) +*---- +* 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(18,18) + DOUBLE PRECISION PIS10,PIS11 +*---- +* 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, 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/ +* + 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(:16)=0.0D0 +*---- +* COMPUTE CORDE =4*V*SIGMA/S=SQRT(3)*COTE*SIGMA (AVERAGE CORDE) +*---- + S2S3=SQRT2*SQRT3 + S3DS2=SQRT3/SQRT2 + CAUX=SQRT3*COTE + CORDE=CAUX*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) + AUX=SQRT3*COSFI + COSA=AUX-SINFI + COSB=AUX+SINFI + AUX=SQRT3*SINFI + SINA=AUX+COSFI + SINB=COSFI-AUX +*---- +* WEIGHTS TIMES BICKLEY NAYLOR FUNCTIONS +*---- + W006=SINFI*COSA + W106=W006*COSB + W005=COSFI*SINA + W105=W005*COSB + W204=COSFI*SINB + W114=SINFI*SINFI*SINB + W224=COSFI*W204 +*---- +* 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(13)=P(13)+WEIGHT*SINB*FKI3(1) + P(14)=P(14)+WEIGHT*W204*FKI4(1) + P(16)=P(16)+WEIGHT*W114*FKI5(1) + P(15)=P(15)+WEIGHT*W224*FKI5(1) + AUX5=WEIGHT*W005 + P(7)=P(7)+AUX5*(FKI4(IAUX1)-FKI4(1)) + P(9)=P(9)+WEIGHT*W105*(FKI5(IAUX1)-FKI5(1)) + AUX=AUX5*(FKI6(IAUX1)-FKI6(1)) + P(11)=P(11)+AUX + P(10)=P(10)+W005*AUX + GO TO 30 +* + 25 WEIGHT=WGAUS(I) + AUX5=WEIGHT*W005 + P(7)=P(7)+AUX5*FKI4(IAUX1) + P(9)=P(9)+WEIGHT*W105*FKI5(IAUX1) + AUX=AUX5*FKI6(IAUX1) + P(11)=P(11)+AUX + P(10)=P(10)+W005*AUX + 30 P(1)=P(1)+WEIGHT*W006*(BI4(0)-FKI4(IAUX2)) + P(3)=P(3)+WEIGHT*W106*(BI5(0)-FKI5(IAUX2)) + AUX=W006*WEIGHT*(.533333333333333-FKI6(IAUX2)) + P(5)=P(5)+AUX + P(4)=P(4)+AUX*W006 + 50 CONTINUE +*---- +* NORMALIZATION +*---- + CORDI=1./CORDE + X1=CORDI/3. + X2=CORDI*SQRT3 + X3=3.*CORDI + X4=X2/SQRT2 + P(1)=X1*P(1) + P(3)=X2*P(3)/6. + P(5)=-X4*P(5) + P(4)=X3*P(4) + P(7)=X1*P(7) + P(9)=X1*P(9)*.5 + P(11)=-X4*P(11) + P(10)=X3*P(10) + P(13)=P(13)/3. + P(14)=SQRT2*P(14) + P(16)=4.*P(16)/3. + P(15)=6.*P(15) + AUX=2.*SQRT2 + P(2)=S3DS2*P(3)-AUX*P(1) + P(8)=3.*S3DS2*P(9)-AUX*P(7) + P(14)=P(14)-AUX*P(13) + COEF=1./(6.*CORDE) + PIS10=(1.-2.*(P(1)+P(7))-P(13))*COEF + PIS(1)=REAL(PIS10) + PIS11=-(2.*(P(2)+P(8))+P(14))*COEF + PIS(2)=REAL(PIS11) + PII=(1.-6.*PIS(1))/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) + P(5)=-S2S3*(1.125*AUX0-.5*TAU0*AUX1+.375*TAU02*AUX2) + P(3)=(-(1.5*ALOGX+1.-SQRT3)*TAU0*.25+(1.-TAU02)/9+ + + AUX*TAU02/3.)*SQRT3 + P(4)=2.25*(1.25*SQRT3-2.)-TAU0*(2.-3.*AUX)+TAU02* + + (2.25*ALOGX-1.5) + XAUX=5.*SQRT3-9. + XAUX0=SQRT3-1.5 + AUX1=AUX*(2.5*ALOG3-4.*ALOG2)-.5 + AUX2=XAUX-2.*ALOG3+.5*ALOGX + P(7)=XAUX0+TAU0*AUX1-TAU02*AUX2/3. + P(11)=-S2S3*(1.125*XAUX0+TAU0*AUX1-.25*AUX2*TAU02) + P(9)=SQRT3/9.+.25*(XAUX-.5*SQRT3*(2.*ALOG3-ALOGX))*TAU0- + + (9.*ALOG3-16.*ALOG2-3.)*TAU02/(3.*PI) + P(10)=2.25*(SQRT3-.75)+TAU0*(3.*AUX-6.)-1.5*TAU02* + + (9.*(2.-SQRT3)+1.5*ALOGX-6.*ALOG3) + AUX0=2.*AUX0 + AUX1=AUX*(.5*ALOG3-ALOG2) + P(13)=AUX0-8.*(1./6.+AUX1)*TAU0-4.*TAU02*(AUX0-.5*ALOG3) + P(15)=4.5*(2.5-SQRT3)-8.*TAU0+36.*TAU02*AUX0 + AUX2=XAUX-2.*ALOG3+.5*ALOGX + P(16)=3.5-2.*SQRT3-64.*TAU0*(AUX1+1./12.)/3.+8.*TAU02* + + (.5*ALOG3-2.*AUX0) + PSQ3=SQRT3/PI + P(14)=SQRT2*(2./3.-6.*TAU0*AUX0+TAU02*(4.+24.*PSQ3* + + (.5*ALOG3-ALOG2))) + PII=COTE*SQRT3*(4.*SQRT3-8.-2.*ALOGX+10.*ALOG3)/12. + PIS(1)=(1.-SIGT*PII)/6. + PIS(2)=-SQRT2*(2.5-2.25*ALOG3+TAU0*(9.+PSQ3*(4.-8.*ALOG2+ + + 3.*ALOG3)-8./SQRT3-(20.*ALOG3-4.*ALOGX)/3.))/12. + AUX=2.*SQRT2 + P(2)=S3DS2*P(3)-AUX*P(1) + P(8)=3.*S3DS2*P(9)-AUX*P(7) + P(14)=P(14)-AUX*P(13) +* + 350 CONTINUE +*---- +* COMPUTE REMAINING PROBABILITIES +*---- + P(4)=P(4)-4.*SQRT3*P(3)+8.*P(1) + P(5)=P(5)+AUX*P(3) + P(6)=(S2S3*P(5)+8.*(-SQRT3*P(3)+P(1))-P(4))*2./9. + P(10)=P(10)-4.*SQRT2*P(8)-8.*P(7) + P(11)=P(11)+AUX*P(9) + P(12)=(8.*(P(7)-SQRT3*P(9))-S2S3*P(11)-P(10))*2./9. + P(15)=P(15)-4.*SQRT2*P(14)-8.*P(13) + P(3)=-P(3) + P(5)=-P(5) + P(9)=-P(9) + P(11)=-P(11) +*---- +* TRANSMISSION MATRIX +*---- + DO 59 I=1,18 + DO 58 J=1,18 + IB=IROT(J,I) + IF(IB.LT.0) THEN + PSS(I,J)=-REAL(P(-IB)) + ELSEIF(IB.GT.0) THEN + PSS(I,J)=REAL(P(IB)) + ELSE + PSS(I,J)=0. + ENDIF + 58 CONTINUE + 59 CONTINUE +*---- +* LEAKAGE PRABABILITIES MATRIX +*---- + PIS(3)=0. + K=3 + DO 56 I=1,5 + K=K+3 + PIS(K-2)=PIS(1) + PIS(K-1)=PIS(2) + PIS(K)=0. + 56 CONTINUE +* + RETURN + END -- cgit v1.2.3