diff options
Diffstat (limited to 'Dragon/src/SYBHTK.f')
| -rw-r--r-- | Dragon/src/SYBHTK.f | 512 |
1 files changed, 512 insertions, 0 deletions
diff --git a/Dragon/src/SYBHTK.f b/Dragon/src/SYBHTK.f new file mode 100644 index 0000000..109966f --- /dev/null +++ b/Dragon/src/SYBHTK.f @@ -0,0 +1,512 @@ +*DECK SYBHTK + SUBROUTINE SYBHTK (NA,NX,NREG,SIDE,RAYRE,ILIGN,INORM,IQW,LR,Z,LI, + 1 IZ,PREC) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute the tracking information related to an hexagonal heterogeneous +* cell. +* +*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 +* NA number of angles in (0,$\\pi$/6). +* NX number of tracks in each sub domain for a given angle. +* NREG number of regions in the cell. +* SIDE side of an hexagon. +* RAYRE radius of each cylinder (RAYRE(1)=0.0). +* ILIGN tracking print flag (=1 to print the tracking). +* INORM track normalization flag (=1 to avoid track normalization). +* IQW equal weight quadrature flag (=1 to use equal weight +* quadratures in angle and space). +* +*Parameters: output +* LR exact size of array Z with +* L.LE.4+3*NA*(13+2*(NREG+1)*NX*NREG). +* Z real tracking information. +* Z(1) to Z(4) contain the numerical orthonormalization +* factors. +* LI size of array IZ with +* L.LE.NREG+4+3*NA*(2+(NREG+1)*(3+2*NREG)). +* IZ integer tracking information. +* IZ(1)=5, IZ(2)=NREG+1 and IZ(3)=3 for an hexagonal cell. +* PREC accuracy obtained if the non-normalized tracks are used +* to integrate the volumes. +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NA,NX,NREG,ILIGN,INORM,IQW,LR,LI,IZ(*) + REAL SIDE,RAYRE(NREG),Z(*),PREC +*---- +* LOCAL VARIABLES +*---- + PARAMETER (PIO2=1.570796327,PI=3.14159265358979,SQ3=1.73205080757) + REAL ZX(64),WX(64),ZA(64),WA(64),ZXJ(64),WXJ(64) + REAL, ALLOCATABLE, DIMENSION(:) :: VAP +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(VAP(NREG)) +* + NA3=3*NA + IF(NX.GT.10) CALL XABORT('SYBHTK: NX IS GREATER THAN 10.') + IF(NA.GT.64) CALL XABORT('SYBHTK: NA IS GREATER THAN 64.') + IF(2.0*RAYRE(NREG).GT.SQ3*SIDE) CALL XABORT('SYBHTK: A RADIUS IS' + 1 //' GREATER THAN HALF A SIDE.') +*---- +* COMPUTE VOLUMES +*---- + VOL=1.5*SQ3*SIDE**2 + DO 10 IR=NREG,1,-1 + R2=PI*RAYRE(IR)**2 + Z(IR)=VOL-R2 + VOL=R2 + 10 CONTINUE +* + IF(IQW.EQ.0) THEN +* GAUSS-LEGENDRE AND GAUSS-JACOBI INTEGRATION POINTS. + CALL ALGPT(NX,-1.,1.,ZX(1),WX(1)) + CALL ALGJP(NX,ZXJ,WXJ) + CALL ALGPT(NA,-1.,-1./3.,ZA(1),WA(1)) + CALL ALGPT(NA,-1./3.,1./3.,ZA(NA+1),WA(NA+1)) + CALL ALGPT(NA,1./3.,1.,ZA(2*NA+1),WA(2*NA+1)) + ELSE +* EQUAL WEIGHT INTEGRATION POINTS. + DO 15 I=1,NX + ZX(I)=(2.0*REAL(I)-1.0)/REAL(NX)-1.0 + WX(I)=2.0/REAL(NX) + ZXJ(I)=0.5*(2.0*REAL(I)-1.0)/REAL(NX) + WXJ(I)=ZXJ(I)/REAL(NX) + 15 CONTINUE + DO 20 I=1,NA3 + ZA(I)=(2.0*REAL(I)-1.0)/REAL(NA3)-1.0 + WA(I)=2.0/REAL(NA3) + 20 CONTINUE + ENDIF + IZ(1)=5 + IZ(2)=NREG+1 + IZ(3)=3 + IZ(4)=NA3 + PREC=0.0 + LI=4 + LR=NREG+4 +*---- +* INTEGRATION IN ANGLE FROM 0 TO PI/2 +*---- + ZN1=0.0 + ZN2=0.0 + ZN3=0.0 + DO 350 IA=1,NA3 + PHI=0.5*PIO2*(ZA(IA)+1.0) + SI=SIN(PHI) + CO=COS(PHI) + TA=SI/CO + FACT1=(SQ3/TA)/(SQ3/TA-1.0) + FACT2=(SQ3/TA)/(SQ3/TA+1.0) + ZN1=ZN1+SI*WA(IA) + ZN2=ZN2+SI*SI*WA(IA) + ZN3=ZN3+SI*SI*SI*WA(IA) + Z(LR+1:LR+6)=0.0 + Z(LR+9)=SI + Z(LR+10)=CO + IF(PHI.LT.PI/6.) THEN + Z(LR+11)=COS(PHI+PI/6.) + Z(LR+12)=SIN(PHI+PI/6.) + ELSE IF(PHI.LT.PI/3.) THEN + Z(LR+1)=COS(PHI-PI/6.) + Z(LR+2)=SIN(PHI-PI/6.) + ELSE + Z(LR+3)=SI + Z(LR+4)=CO + ENDIF + Z(LR+13)=WA(IA) + LR=LR+13 +*---- +* FIRST ANGULAR DOMAIN (0 TO PI/6) +*---- + L4=LI+1 + IZ(LI+1)=0 + IZ(LI+2)=0 + LI=LI+2 + IF(PHI.GT.PIO2/3.0) GO TO 120 + X1=0.0 + XLIM=MIN(SIDE,0.5*SIDE*(SQ3/TA-1.0)) + DLIM=0.5*SIDE*SQ3*CO+(0.5*SIDE-XLIM)*SI + DO 100 K0=NREG,1,-1 + KMAX=NREG-K0+1 + X2=MIN(XLIM,XLIM-(RAYRE(K0)-DLIM)/SI) + L3=LR+1 + L5=LI+1 + LI=LI+3 + VAP(:NREG)=0.0 + DO 50 IX=1,NX + IF(K0.EQ.NREG) THEN + S=0.5*(X2-X1)*SI*WX(IX) + X=X1+0.5*(X2-X1)*(1.0+ZX(IX)) + ELSE +* FLURIG CHANGE OF VARIABLE. + S=2.0*(X2-X1)*SI*WXJ(IX) + X=X1+(X2-X1)*ZXJ(IX)**2 + ENDIF + Z(LR+1)=S*WA(IA) + LR=LR+1 + C=0.5*SIDE*SQ3*SI-(0.5*SIDE-X)*CO + D=0.5*SIDE*SQ3*CO+(0.5*SIDE-X)*SI + D=D*D + SUM=0.0 + CORDE=0.0 + DO 30 K=NREG,K0+1,-1 + RR=RAYRE(K)**2-D + CORDE=SQRT(RR) + DEL=C-CORDE + SUM=SUM+DEL + Z(LR+NREG-K+1)=DEL + VAP(K)=VAP(K)+DEL*S + C=CORDE + 30 CONTINUE + IF(KMAX.NE.1) THEN + DEL=2.0*CORDE + SUM=SUM+DEL + Z(LR+KMAX)=DEL + VAP(K)=VAP(K)+DEL*S + DO 40 I=1,KMAX-2 + DEL=Z(LR+KMAX-I) + SUM=SUM+DEL + Z(LR+KMAX+I)=DEL + VAP(K+I)=VAP(K+I)+DEL*S + 40 CONTINUE + ENDIF + LR=LR+2*KMAX-1 + DEL=X*FACT1/CO-SUM + Z(LR)=DEL + VAP(NREG)=VAP(NREG)+DEL*S + 50 CONTINUE + DO 60 K=KMAX-1,1-KMAX,-1 + IZ(LI+K+KMAX)=5+ABS(K)+1+NREG-KMAX + 60 CONTINUE + LI=LI+2*KMAX + IZ(L5)=2*KMAX-1 + IZ(L5+1)=NX + IZ(L5+2)=4 ! ISURF + IZ(LI)=5 ! JSURF +*---- +* VOLUME NORMALIZATION +*---- + IF((INORM.EQ.0).AND.(K0.LT.NREG)) THEN + DLIM1=0.5*SIDE*SQ3*CO+(0.5*SIDE-X2)*SI + DLIM2=0.5*SIDE*SQ3*CO+(0.5*SIDE-X1)*SI + VW1=0.0 + SUMVAP=0.0 + DO 70 I=K0,NREG-1 + SUMVAP=SUMVAP+VAP(I) + RW=RAYRE(I+1) + VEX1=RW*RW*ACOS(DLIM1/RW)-DLIM1*SQRT(RW*RW-DLIM1*DLIM1) + IF(RW.GT.DLIM2) + 1 VEX1=VEX1-(RW*RW*ACOS(DLIM2/RW)-DLIM2*SQRT(RW*RW-DLIM2*DLIM2)) + SUM=(VEX1-VW1)/VAP(I) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX1-VW1)/(1.5*SQ3*SIDE**2)) + VW1=VEX1 + VAP(I)=SUM + 70 CONTINUE + VEX1=0.5*(SIDE*SQ3*SI-(SIDE-X1-X2)*CO)*(X2-X1)*SI + VEX2=0.5*FACT1*TA*(X2*X2-X1*X1)-VEX1 + SUM=(VEX1-0.5*VW1)/(VEX1-0.5*SUMVAP) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX1-0.5*VW1)/(1.5*SQ3*SIDE**2)) + VEX1=SUM + SUM=(VEX2-0.5*VW1)/(VEX2-0.5*SUMVAP) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX2-0.5*VW1)/(1.5*SQ3*SIDE**2)) + VEX2=SUM + DO 90 IX=1,NX + KMAX=(IZ(L5)+1)/2 + Z(L3+KMAX)=Z(L3+KMAX)*VAP(K0) + DO 80 I=1,KMAX-2 + Z(L3+KMAX-I)=Z(L3+KMAX-I)*VAP(K0+I) + Z(L3+KMAX+I)=Z(L3+KMAX+I)*VAP(K0+I) + 80 CONTINUE + Z(L3+1)=Z(L3+1)*VEX1 + Z(L3+2*KMAX-1)=Z(L3+2*KMAX-1)*VEX2 + L3=L3+2*KMAX + 90 CONTINUE + ENDIF + IZ(L4)=IZ(L4)+1 + IF(X2.GE.XLIM) GO TO 120 + X1=X2 + 100 CONTINUE +*---- +* SECOND ANGULAR DOMAIN (PI/6 TO PI/3) +*---- + 120 IF((PHI.LE.PI/6.0).OR.(PHI.GT.2.0*PIO2/3.0)) GO TO 240 + X1=0.5*SIDE*(SQ3/TA-1.0) + XLIM=SIDE + DLIM=0.5*SIDE*SQ3*CO+(0.5*SIDE-XLIM)*SI + DO 230 K0=NREG,1,-1 + KMAX=NREG-K0+1 + X2=MIN(XLIM,XLIM-(RAYRE(K0)-DLIM)/SI) + IF(X2.LE.X1) GO TO 230 + L3=LR+1 + L5=LI+1 + LI=LI+3 + VAP(:NREG)=0.0 + DO 150 IX=1,NX + IF(K0.EQ.NREG) THEN + S=0.5*(X2-X1)*SI*WX(IX) + X=X1+0.5*(X2-X1)*(1.0+ZX(IX)) + ELSE +* FLURIG CHANGE OF VARIABLE. + S=2.0*(X2-X1)*SI*WXJ(IX) + X=X1+(X2-X1)*ZXJ(IX)**2 + ENDIF + Z(LR+1)=S*WA(IA) + LR=LR+1 + C=0.5*SIDE*SQ3*SI-(0.5*SIDE-X)*CO + D=0.5*SIDE*SQ3*CO+(0.5*SIDE-X)*SI + D=D*D + SUM=0.0 + CORDE=0.0 + DO 130 K=NREG,K0+1,-1 + RR=RAYRE(K)**2-D + CORDE=SQRT(RR) + DEL=C-CORDE + SUM=SUM+DEL + Z(LR+NREG-K+1)=DEL + VAP(K)=VAP(K)+DEL*S + C=CORDE + 130 CONTINUE + IF(KMAX.NE.1) THEN + DEL=2.0*CORDE + SUM=SUM+DEL + Z(LR+KMAX)=DEL + VAP(K)=VAP(K)+DEL*S + DO 140 I=1,KMAX-2 + DEL=Z(LR+KMAX-I) + SUM=SUM+DEL + Z(LR+KMAX+I)=DEL + VAP(K+I)=VAP(K+I)+DEL*S + 140 CONTINUE + ENDIF + LR=LR+2*KMAX-1 + DEL=(X+SIDE)*FACT2/CO-SUM + Z(LR)=DEL + VAP(NREG)=VAP(NREG)+DEL*S + 150 CONTINUE + DO 160 K=KMAX-1,1-KMAX,-1 + IZ(LI+K+KMAX)=5+ABS(K)+1+NREG-KMAX + 160 CONTINUE + LI=LI+2*KMAX + IZ(L5)=2*KMAX-1 + IZ(L5+1)=NX + IZ(L5+2)=4 ! ISURF + IZ(LI)=0 ! JSURF +*---- +* VOLUME NORMALIZATION +*---- + IF((INORM.EQ.0).AND.(K0.LT.NREG)) THEN + DLIM1=0.5*SIDE*SQ3*CO+(0.5*SIDE-X2)*SI + DLIM2=0.5*SIDE*SQ3*CO+(0.5*SIDE-X1)*SI + VW1=0.0 + SUMVAP=0.0 + DO 200 I=K0,NREG-1 + SUMVAP=SUMVAP+VAP(I) + RW=RAYRE(I+1) + VEX1=RW*RW*ACOS(DLIM1/RW)-DLIM1*SQRT(RW*RW-DLIM1*DLIM1) + IF(RW.GT.DLIM2) + 1 VEX1=VEX1-(RW*RW*ACOS(DLIM2/RW)-DLIM2*SQRT(RW*RW-DLIM2*DLIM2)) + SUM=(VEX1-VW1)/VAP(I) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX1-VW1)/(1.5*SQ3*SIDE**2)) + VW1=VEX1 + VAP(I)=SUM + 200 CONTINUE + VEX1=0.5*(SIDE*SQ3*SI-(SIDE-X1-X2)*CO)*(X2-X1)*SI + VEX2=0.5*FACT2*TA*(X2-X1)*(X1+X2+2.0*SIDE)-VEX1 + SUM=(VEX1-0.5*VW1)/(VEX1-0.5*SUMVAP) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX1-0.5*VW1)/(1.5*SQ3*SIDE**2)) + VEX1=SUM + SUM=(VEX2-0.5*VW1)/(VEX2-0.5*SUMVAP) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX2-0.5*VW1)/(1.5*SQ3*SIDE**2)) + VEX2=SUM + DO 220 IX=1,NX + KMAX=(IZ(L5)+1)/2 + Z(L3+KMAX)=Z(L3+KMAX)*VAP(K0) + DO 210 I=1,KMAX-2 + Z(L3+KMAX-I)=Z(L3+KMAX-I)*VAP(K0+I) + Z(L3+KMAX+I)=Z(L3+KMAX+I)*VAP(K0+I) + 210 CONTINUE + Z(L3+1)=Z(L3+1)*VEX1 + Z(L3+2*KMAX-1)=Z(L3+2*KMAX-1)*VEX2 + L3=L3+2*KMAX + 220 CONTINUE + ENDIF + IZ(L4)=IZ(L4)+1 + X1=X2 + 230 CONTINUE +*---- +* THIRD ANGULAR DOMAIN (PI/3 TO PI/2) +*---- + 240 IF(PHI.LE.2.0*PIO2/3.0) GO TO 350 + X1=SQ3*SIDE/TA + XLIM=0.5*(SIDE+X1) + DO 340 K0=NREG,1,-1 + KMAX=NREG-K0+1 + X2=XLIM-RAYRE(K0)/SI + IF(X2.LE.X1) GO TO 340 + L3=LR+1 + L5=LI+1 + LI=LI+3 + VAP(:NREG)=0.0 + DO 270 IX=1,NX + IF(K0.EQ.NREG) THEN + S=0.5*(X2-X1)*SI*WX(IX) + X=X1+0.5*(X2-X1)*(1.0+ZX(IX)) + ELSE +* FLURIG CHANGE OF VARIABLE. + S=2.0*(X2-X1)*SI*WXJ(IX) + X=X1+(X2-X1)*ZXJ(IX)**2 + ENDIF + Z(LR+1)=S*WA(IA) + LR=LR+1 + C=0.5*SIDE*SQ3*SI-(0.5*SIDE-X)*CO + D=0.5*SIDE*SQ3*CO+(0.5*SIDE-X)*SI + D=D*D + SUM=0.0 + CORDE=0.0 + DO 250 K=NREG,K0+1,-1 + RR=RAYRE(K)**2-D + CORDE=SQRT(RR) + DEL=C-CORDE + SUM=SUM+DEL + Z(LR+NREG-K+1)=DEL + VAP(K)=VAP(K)+DEL*S + C=CORDE + 250 CONTINUE + IF(KMAX.NE.1) THEN + DEL=2.0*CORDE + SUM=SUM+DEL + Z(LR+KMAX)=DEL + VAP(K)=VAP(K)+DEL*S + DO 260 I=1,KMAX-2 + DEL=Z(LR+KMAX-I) + SUM=SUM+DEL + Z(LR+KMAX+I)=DEL + VAP(K+I)=VAP(K+I)+DEL*S + 260 CONTINUE + ENDIF + IF(KMAX.NE.KMAX) CALL XABORT('BUG') + LR=LR+2*KMAX-1 + DEL=SQ3*SIDE/SI-SUM + Z(LR)=DEL + VAP(NREG)=VAP(NREG)+DEL*S + 270 CONTINUE + DO 280 K=KMAX-1,1-KMAX,-1 + IZ(LI+K+KMAX)=5+ABS(K)+1+NREG-KMAX + 280 CONTINUE + LI=LI+2*KMAX + IZ(L5)=2*KMAX-1 + IZ(L5+1)=NX + IZ(L5+2)=4 ! ISURF + IZ(LI)=1 ! JSURF +*---- +* VOLUME NORMALIZATION +*---- + IF((INORM.EQ.0).AND.(K0.LT.NREG)) THEN + DLIM1=0.5*SIDE*SQ3*CO+(0.5*SIDE-X2)*SI + DLIM2=0.5*SIDE*SQ3*CO+(0.5*SIDE-X1)*SI + VW1=0.0 + SUMVAP=0.0 + DO 310 I=K0,NREG-1 + SUMVAP=SUMVAP+VAP(I) + RW=RAYRE(I+1) + VEX1=RW*RW*ACOS(DLIM1/RW)-DLIM1*SQRT(RW*RW-DLIM1*DLIM1) + IF(RW.GT.DLIM2) + 1 VEX1=VEX1-(RW*RW*ACOS(DLIM2/RW)-DLIM2*SQRT(RW*RW-DLIM2*DLIM2)) + SUM=(VEX1-VW1)/VAP(I) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX1-VW1)/(1.5*SQ3*SIDE**2)) + VW1=VEX1 + VAP(I)=SUM + 310 CONTINUE + VEX1=0.5*(SIDE*SQ3*SI-(SIDE-X1-X2)*CO)*(X2-X1)*SI + VEX2=(X2-X1)*SQ3*SIDE-VEX1 + SUM=(VEX1-0.5*VW1)/(VEX1-0.5*SUMVAP) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX1-0.5*VW1)/(1.5*SQ3*SIDE**2)) + VEX1=SUM + SUM=(VEX2-0.5*VW1)/(VEX2-0.5*SUMVAP) + PREC=MAX(PREC,ABS(1.0-SUM)*(VEX2-0.5*VW1)/(1.5*SQ3*SIDE**2)) + VEX2=SUM + DO 330 IX=1,NX + KMAX=(IZ(L5)+1)/2 + Z(L3+KMAX)=Z(L3+KMAX)*VAP(K0) + DO 320 I=1,KMAX-2 + Z(L3+KMAX-I)=Z(L3+KMAX-I)*VAP(K0+I) + Z(L3+KMAX+I)=Z(L3+KMAX+I)*VAP(K0+I) + 320 CONTINUE + Z(L3+1)=Z(L3+1)*VEX1 + Z(L3+2*KMAX-1)=Z(L3+2*KMAX-1)*VEX2 + L3=L3+2*KMAX + 330 CONTINUE + ENDIF + IZ(L4)=IZ(L4)+1 + X1=X2 + 340 CONTINUE + 350 CONTINUE + ZN1=0.5*ZN1*PIO2 + ZN2=0.5*ZN2*PIO2 + ZN3=0.5*ZN3*PIO2 + Z(NREG+1)=1.0/SQRT(ZN1) + Z(NREG+2)=1.0/SQRT(0.75*ZN3-0.7205061948*ZN2*ZN2/ZN1) + Z(NREG+3)=Z(NREG+2)*0.8488263632*ZN2/ZN1 + Z(NREG+4)=2.0/SQRT(3.0*(ZN1-ZN3)) +*---- +* TRACKING INFORMATION OUTPUT +*---- + IF(ILIGN.EQ.1) THEN + L1I=IZ(1)-1 + L1R=IZ(2)-1 + WRITE(6,500) (Z(L1R+I),I=1,4) + L1R=L1R+4 + L2=0 + DO 375 IA=1,NA3 + MNT=IZ(L1I+1) + L1I=L1I+2 + ZSIN=Z(L1R+9) + ZCOS=Z(L1R+10) + L1R=L1R+13 + DO 370 IMNT=1,MNT + NH=IZ(L1I+1) + NX=IZ(L1I+2) + L1I=L1I+3 + DO 360 IX=1,NX + L2=L2+1 + IF((IMNT.EQ.1).AND.(IX.EQ.1)) THEN + WRITE(6,510) L2,ZSIN,ZCOS,Z(L1R+1),NH,(Z(L1R+I+1),I=1,NH) + ELSE + WRITE(6,520) L2,Z(L1R+1),NH,(Z(L1R+I+1),I=1,NH) + ENDIF + L1R=L1R+NH+1 + 360 CONTINUE + L1I=L1I+NH+1 + 370 CONTINUE + 375 CONTINUE + ENDIF +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(VAP) + RETURN +* + 500 FORMAT (1H1//30H TRACKING INFORMATION LISTING.//12H NUMERICAL O, + 1 27HRTHONORMALIZATION FACTORS =,1P,4E12.4//6H TRACK) + 510 FORMAT (1X,I5,7H SIN =,1P,E10.3,7H COS =,E10.3,9H WEIGHT =, + 1 E10.3,6H NH =,I3,12H SEGMENTS =,5E10.3:/(80X,5E10.3)) + 520 FORMAT (1X,I5,34X,9H WEIGHT =,1P,E10.3,6H NH =,I3,10H SEGMENTS, + 1 2H =,5E10.3:/(80X,5E10.3)) + END |