summaryrefslogtreecommitdiff
path: root/Dragon/src/SYB7T0.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SYB7T0.f')
-rw-r--r--Dragon/src/SYB7T0.f387
1 files changed, 387 insertions, 0 deletions
diff --git a/Dragon/src/SYB7T0.f b/Dragon/src/SYB7T0.f
new file mode 100644
index 0000000..007c983
--- /dev/null
+++ b/Dragon/src/SYB7T0.f
@@ -0,0 +1,387 @@
+*DECK SYB7T0
+ SUBROUTINE SYB7T0(MNA,NRD,COTE,RAYONS,JMINR,XCOTE,LFAIRE,DELR,
+ 1 IQW,PWA2,ZWA2,NXMIN,NXMAX,NZR,ZZR,NZI,ZZI)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the tracking information related to an hexagonal sectorized
+* heterogeneous cell.
+*
+*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
+* MNA number of angles in (0,$\\pi$/6).
+* NRD one plus the number of tubes in the cell.
+* COTE length of the hexagon side.
+* RAYONS radius of each cylinder.
+* JMINR first interception with side.
+* XCOTE interceptions with side.
+* LFAIRE tracking calculation flag (=.FALSE. only compute the number
+* of real and integer tracking elements).
+* DELR half distance between the tracks.
+* IQW equal weight quadrature flag (=1 to use equal weight
+* quadratures in angle and space).
+* PWA2 weights of the angular quadrature set.
+* ZWA2 base points of the angular quadrature set.
+*
+*Parameters: output
+* NXMIN minimum number of tracks per region.
+* NXMAX maximum number of tracks per region.
+* NZR number of real tracking elements.
+* ZZR real tracking information.
+* NZI number of integer tracking elements.
+* ZZI integer tracking information.
+*
+*-----------------------------------------------------------------------
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER MNA,NRD,JMINR,IQW,NXMIN,NXMAX,NZR,NZI,ZZI(*)
+ REAL COTE,RAYONS(NRD-1),XCOTE(NRD),DELR,PWA2(64),
+ & ZWA2(64),ZZR(*)
+ LOGICAL LFAIRE
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (PI314 = 3.141592653589793)
+ PARAMETER (PI6 = PI314 / 6)
+ PARAMETER (PI12 = PI314 / 12)
+ PARAMETER (SQRT3 = 1.732050807568877)
+ PARAMETER (SQRT32 = SQRT3 / 2)
+ REAL ANGLES(5)
+ REAL COSECT(3)
+ LOGICAL LGTRAE
+ LOGICAL LGTRAW
+ LOGICAL LVERIF
+ LOGICAL NTRIOR
+ REAL ORIGIN(4)
+ REAL PENTES(4)
+ CHARACTER * 4 TYSUIT
+ REAL WX(64)
+ REAL ZX(64)
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: IXRAYO
+ REAL, ALLOCATABLE, DIMENSION(:) :: HXRAYO
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(IXRAYO(NRD*2+3),HXRAYO(NRD*2+3))
+*----
+* Interpolation des Trajectoires
+*----
+ IZI = 0
+ IZR = 0
+ LVERIF = .FALSE.
+ DELR2 = DELR * 2.
+ HAUTEU = COTE * SQRT32
+*----
+* /Debut/ Boucle sur les Angles
+*----
+ DO 350 IA = 1, MNA
+ MNT = 0
+ IZI = IZI + 1
+ IZIT = IZI
+*
+ WANGIA = PWA2(IA)
+ ZAIA = ZWA2(IA) + 1.
+ YAIA = 1. - ZWA2(IA)
+ WANGIA = WANGIA / 6.
+*
+ PHI = PI12 * ZAIA
+ PHI6 = PI12 * YAIA
+ PHI3 = PHI + PI6
+ COSPHI = COS(PHI)
+ SINPHI = SIN(PHI)
+ COSPH3 = COS(PHI3)
+ SINPH3 = SIN(PHI3)
+ COSPH6 = COS(PHI6)
+ SINPH6 = SIN(PHI6)
+ TANPH6 = TAN(PHI6)
+ IF(LFAIRE) THEN
+ ZZR(IZR+1) = COSPH6
+ ZZR(IZR+2) = SINPH6
+ ZZR(IZR+3) = SINPHI
+ ZZR(IZR+4) = COSPHI
+ ZZR(IZR+5) = COSPH3
+ ZZR(IZR+6) = -SINPH3
+ ZZR(IZR+7) = WANGIA
+ ENDIF
+ IZR = IZR + 7
+ ANGLES(1) = COSPHI
+ ANGLES(2) = SINPHI
+ ANGLES(3) = - TAN(PHI3)
+ ANGLES(4) = TANPH6
+ ANGLES(5) = COSPHI / SINPHI
+*
+ PENTES(1) = TANPH6
+ PENTES(2) = ANGLES(5)
+ PENTES(3) = - TAN(PHI3)
+ PENTES(4) = PENTES(1)
+ ORIGIN(1) = - HAUTEU / COS(PHI6)
+ ORIGIN(2) = - HAUTEU / SINPHI
+ ORIGIN(3) = HAUTEU / COS(PHI3)
+ ORIGIN(4) = - ORIGIN(1)
+*
+ COSECT(1) = COS(PHI3)
+ COSECT(2) = COS(PHI6)
+ COSECT(3) = SINPHI
+*----
+* Rayon Max Pour cet angle(Rext)
+*----
+ XIZERO = HAUTEU * TANPH6
+ DO MRA = JMINR, NRD-1
+ IF(XCOTE(MRA) .GE. XIZERO) GOTO 140
+ ENDDO
+ MRA = NRD
+ 140 MRAE = MRA
+ MRAW = MRA
+*----
+* Les Lignes d'integrations sont limites
+* et par les secteurs
+* et par les couronnes
+*
+* Demarrage au centre de l'hexagone
+*----
+ IHMAX = NRD - MRA
+ IHMIN = IHMAX + 1
+ HXRAYO(IHMIN) = ORIGIN(1)
+ DO IR = MRAW, 2, -1
+ IHMAX = IHMAX + 1
+ IXRAYO(IHMAX) = IR
+ HXRAYO(IHMAX+1) = - RAYONS(IR-1)
+ ENDDO
+ DO IS = 1, 3
+ IXRAYO(IHMAX+IS) = 1
+ HXRAYO(IHMAX+IS+1) = 0.
+ ENDDO
+ IHMAX = IHMAX + 3
+ HXRAYO(IHMAX+1) = 0.
+ DO IR = 1, MRAE-1
+ IHMAX = IHMAX + 1
+ IXRAYO(IHMAX) = IR
+ HXRAYO(IHMAX+1) = RAYONS(IR)
+ ENDDO
+ IHMAX = IHMAX + 1
+ IXRAYO(IHMAX) = MRAE
+ HXRAYO(IHMAX+1) = ORIGIN(4)
+*
+ DELTCW = COTE * COSECT(1)
+ DELTAS = COTE * COSECT(2)
+ DELTCE = COTE * COSECT(3)
+ ORIPHI = HAUTEU * COSPHI
+ ORIPH6 = HAUTEU * SINPH6
+ ORIPH3 = HAUTEU * SINPH3
+*
+ ISW2 = 2
+ CALL SYB7TW(NRD, JMINR, XCOTE, IXRAYO(IHMIN),
+ & DELTCW, DELTAS, ORIPH6, ORIPHI,
+ & COSPH6, SINPHI,
+ & ISW2, LGTRAW, DELTAW)
+ ISW = ISW2 / 2
+*
+ ISE2 = 8
+ CALL SYB7TE(NRD, JMINR, XCOTE, IXRAYO(IHMAX),
+ & DELTCE, DELTAS, ORIPH6, ORIPH3,
+ & COSPH6, COSPH3,
+ & ISE2, LGTRAE, DELTAE)
+*
+ DELTAC = 0.
+ NTRIOR = .TRUE.
+ IRSUIT = 0
+*----
+* /Debut/ Boucle sur les Trajectoires
+*----
+ DO WHILE(NTRIOR)
+ NHMAX = IHMAX + 1 - IHMIN
+ DELTAP = DELTAC
+*
+ TYSUIT = 'Sud'
+ IHSUIT = IHMIN
+ DELTAC = DELTAS
+ IF(DELTAW .LT. DELTAC) THEN
+ TYSUIT = 'West'
+ IHSUIT = IHMIN
+ DELTAC = DELTAW
+ ENDIF
+ IF(DELTAE .LT. DELTAC) THEN
+ TYSUIT = 'Est'
+ IHSUIT = IHMAX
+ DELTAC = DELTAE
+ ENDIF
+ CALL SYB7TN(IHMIN, IHMAX, IXRAYO, ISW,
+ & COSECT, NRD-1, RAYONS,
+ & TYSUIT, IHSUIT, DELTAC, IRSUIT)
+*
+ DELTAX = DELTAC - DELTAP
+ IF(DELTAX .LE. DELR2) THEN
+ NX = 1
+ ZX(1) = 0.0
+ WX(1) = 2.0
+ ELSE
+ NX = INT(DELTAX / DELR2 + 1)
+ IF(IQW.EQ.0) THEN
+* GAUSS-LEGENDRE INTEGRATION POINTS. ZX(I) IS NOT USED.
+ IF(NX.GT.20) THEN
+ IF(NX.LT.24) THEN
+ NX=24
+ ELSE IF(NX.LT.28) THEN
+ NX=28
+ ELSE IF(NX.LT.32) THEN
+ NX=32
+ ELSE IF(NX.LT.64) THEN
+ NX=64
+ ELSE IF(NX.GT.64) THEN
+ CALL XABORT('SYB7T0: GAUSS OVERFLOW.')
+ ENDIF
+ ENDIF
+ CALL ALGPT(NX,-1.0,1.0,ZX,WX)
+ ELSE
+* EQUAL WEIGHT INTEGRATION POINTS.
+ DO 30 I=1,NX
+ ZX(I)=(2.0*REAL(I)-1.0)/REAL(NX)-1.0
+ WX(I)=2.0/REAL(NX)
+ 30 CONTINUE
+ ENDIF
+ ENDIF
+ NXMIN = MIN(NX, NXMIN)
+ NXMAX = MAX(NX, NXMAX)
+*
+ MNT = MNT + 1
+ IF(LFAIRE) THEN
+ ZZI(IZI+1) = NHMAX
+ ZZI(IZI+2) = NX
+ DO I=0,NHMAX-1
+ ZZI(IZI+3+I)=IXRAYO(IHMIN+I)
+ ENDDO
+*
+ DELTAR = DELTAP
+ DO 250 IX = 1, NX
+*
+ WW = 0.5 * WX(IX)
+*
+ IZR = IZR + 1
+ ZZR(IZR) = WW * WANGIA * DELTAX
+ DDELTA = DELTAX * WW
+ DELTAR = DELTAR + DDELTA
+*
+* Position de L'intersection Gauche(West)
+ ZZW = DELTAR * PENTES(ISW) + ORIGIN(ISW)
+*
+* Position de L'intersection Droite(Est)
+ ISE = ISE2 / 2
+ ZZE = DELTAR * PENTES(ISE) + ORIGIN(ISE)
+*----
+* Longueur des intersections
+*----
+ CALL SYB7TC(DELTAR, DDELTA, ANGLES(ISW+2), NHMAX, ZZI(IZI+3),
+ & NRD-1, RAYONS, ZZW, ZZE,
+ & ZZR(IZR+1), HXRAYO(IHMIN))
+ IZR = IZR + NHMAX
+*
+ 250 CONTINUE
+ ELSE
+*----
+* Pour le comptage il faudra compter
+* pour plus tard (les douzes symetries)
+* deux indices de plus pour les surfaces Entrante et Sortante
+*----
+ IZR = IZR +(NHMAX + 1) * NX
+ IZI = IZI + 2
+ ENDIF
+*
+ IZI = IZI + NHMAX + 2
+*
+ IF(TYSUIT .EQ. 'Est') THEN
+ IRC = IXRAYO(IHMAX)
+ IF(LGTRAE) THEN
+ IHMAX = IHMAX - 1
+ IF(IRC .EQ. IXRAYO(IHMAX)) ISE2 = ISE2 - 1
+ ELSEIF(IRC .GE. NRD) THEN
+ NTRIOR = .FALSE.
+ ELSE
+ IHMAX = IHMAX + 1
+ IXRAYO(IHMAX) = IRC + 1
+ HXRAYO(IHMAX+1) = HXRAYO(IHMAX)
+ ENDIF
+ IF(NTRIOR) THEN
+ CALL SYB7TE(NRD, JMINR, XCOTE, IXRAYO(IHMAX),
+ & DELTCE, DELTAS, ORIPH6, ORIPH3,
+ & COSPH6, COSPH3,
+ & ISE2, LGTRAE, DELTAE)
+ ENDIF
+ ELSEIF(TYSUIT .EQ. 'West') THEN
+ IRC = IXRAYO(IHMIN)
+ IF(LGTRAW) THEN
+ IHMIN = IHMIN + 1
+ IF(IRC .EQ. IXRAYO(IHMIN)) ISW2 = ISW2 + 1
+ ELSEIF(IRC .GE. NRD) THEN
+ NTRIOR = .FALSE.
+ ELSE
+ IHMIN = IHMIN - 1
+ IXRAYO(IHMIN) = IRC + 1
+ HXRAYO(IHMIN) = HXRAYO(IHMIN+1)
+ ENDIF
+ NTRIOR = ISW .LE. 2
+ IF(NTRIOR) THEN
+ CALL SYB7TW(NRD, JMINR, XCOTE, IXRAYO(IHMIN),
+ & DELTCW, DELTAS, ORIPH6, ORIPHI,
+ & COSPH6, SINPHI,
+ & ISW2, LGTRAW, DELTAW)
+ IF((ISW2 / 2) .NE. ISW) THEN
+ IF(LFAIRE) THEN
+ ZZI(IZI+1) = 0
+ ZZI(IZI+2) = 0
+ MNT = MNT + 1
+ ENDIF
+ IZI = IZI+2
+ ENDIF
+ ISW = ISW2 / 2
+ IF(ISW2 .EQ. 5) THEN
+ IF(LGTRAW) THEN
+ HC = DELTAW * PENTES(ISW) + ORIGIN(ISW)
+ LGTRAW = HC .GE. 0.
+ ENDIF
+ ENDIF
+ ENDIF
+ ELSEIF(TYSUIT .EQ. 'Coin') THEN
+ IXRAYO(IHSUIT) = IRSUIT
+ ELSEIF(TYSUIT .EQ. 'Tang') THEN
+*
+** Decalage du tableau entier de 2 Cases
+ DO I = 2, IHMAX - IHSUIT
+ IXRAYO(IHSUIT+I-2) = IXRAYO(IHSUIT+I)
+ ENDDO
+*
+** Decalage du tableau reel de 2 Cases
+ DO I = 2, IHMAX + 1 - IHSUIT
+ HXRAYO(IHSUIT+I-2) = HXRAYO(IHSUIT+I)
+ ENDDO
+*
+ IHMAX = IHMAX - 2
+ ELSE
+ NTRIOR = .FALSE.
+ ENDIF
+*
+* /Fin/ Boucle sur les Trajectoires
+ ENDDO
+*
+* /Fin/ Boucle sur les Angles
+ IF(LFAIRE) ZZI(IZIT) = MNT
+ 350 CONTINUE
+ NZI = IZI
+ NZR = IZR
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(HXRAYO,IXRAYO)
+*
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 17:29:58 +0000