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/XHX2D0.f | 207 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 207 insertions(+) create mode 100644 Dragon/src/XHX2D0.f (limited to 'Dragon/src/XHX2D0.f') 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 -- cgit v1.2.3