summaryrefslogtreecommitdiff
path: root/Dragon/src/SYB4TH.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SYB4TH.f')
-rw-r--r--Dragon/src/SYB4TH.f375
1 files changed, 375 insertions, 0 deletions
diff --git a/Dragon/src/SYB4TH.f b/Dragon/src/SYB4TH.f
new file mode 100644
index 0000000..bfb5c97
--- /dev/null
+++ b/Dragon/src/SYB4TH.f
@@ -0,0 +1,375 @@
+*DECK SYB4TH
+ SUBROUTINE SYB4TH (NRMAX,RAYONS,PENTES,ORIGIN,ISCMIN,ISCMAX,
+ & ANGLES,DXMIN,DELTAC,ISCW,ISCE,ISXW,ISXE,NHMAX,IXRAYO,HXRAYO)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the intersection lenghts of a track in a sectorized Cartesian
+* 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
+* NRMAX number of radius.
+* RAYONS radius of each tube.
+* PENTES slope of the track with respect to the axis or to the sides.
+* ORIGIN origin of the track.
+* ISCMIN first sector with a possible intersection.
+* ISCMAX last sector with a possible intersection.
+* ANGLES slope of the sectors.
+* DXMIN accuracy.
+* DELTAC position of the track.
+*
+*Parameters: output
+* ISCW index of the first sector (west/left).
+* ISCE index of the last sector (east/right).
+* ISXW index of the first side (west/left). Equal to 1 or 2.
+* ISXE index of the last side (east/right). Equal to 3 or 4.
+* NHMAX number of intersections.
+* IXRAYO indices of the intersecting tubes.
+* HXRAYO position of the intersections (boundary).
+*
+*-----------------------------------------------------------------------
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER NRMAX,ISCMIN,ISCMAX,ISCW,ISCE,ISXW,ISXE,NHMAX,
+ & IXRAYO(NHMAX)
+ REAL RAYONS(NRMAX),PENTES(4),ORIGIN(4),ANGLES(ISCMAX),DXMIN,
+ & DELTAC,HXRAYO(NHMAX+1)
+*----
+* LOCAL VARIABLES
+*----
+ CHARACTER * 4 AJOU(2)
+ LOGICAL AJRAYO
+ LOGICAL AJSECT
+*
+* Remarque Importante :
+* Pour choisir entre deux lignes Proches
+* Il faut prendre le cas correspondant a
+* DELTAC + DXMIN
+* ZHC + DXMIN
+*
+ DXMAX = 2.
+ D2 = DELTAC * DELTAC
+ DR = DELTAC + DXMIN
+ ZH1 = DR * PENTES(1) + ORIGIN(1)
+ ZH2 = DR * PENTES(2) + ORIGIN(2)
+ IF (ZH2 .GE. (ZH1 - DXMIN)) THEN
+ ZHW = DELTAC * PENTES(2) + ORIGIN(2)
+ ISXW = 2
+ ELSE
+ ZHW = DELTAC * PENTES(1) + ORIGIN(1)
+ ISXW = 1
+ ENDIF
+ ZH3 = DELTAC * PENTES(3) + ORIGIN(3)
+ ZH4 = DELTAC * PENTES(4) + ORIGIN(4)
+ IF (ZH3 .LE. (ZH4 + DXMIN)) THEN
+ ZHE = ZH3
+ ISXE = 3
+ ELSE
+ ZHE = ZH4
+ ISXE = 4
+ ENDIF
+*
+* Recherche du Premier Secteur
+* L'Angle ISC delimite les Secteurs ISC,ISC+1
+ DO 101 ISCW = ISCMIN, ISCMAX
+ ZHS = DELTAC * ANGLES(ISCW)
+ IF (ZHS .GT. (ZHW + DXMIN)) THEN
+ GOTO 102
+ ENDIF
+ 101 CONTINUE
+ ISCW = ISCMAX + 1
+ ZHS = ZHE + DXMAX
+ 102 CONTINUE
+*
+* Couronne Interne
+ DO 201 NRMIN = 1, NRMAX
+ IF (RAYONS(NRMIN) .GT. DR) THEN
+ R2 = RAYONS(NRMIN) * RAYONS(NRMIN)
+ H2 = R2 - D2
+ IF (H2 .GT. 0.) THEN
+ H0 = SQRT(H2)
+ GOTO 202
+ ENDIF
+ ENDIF
+ 201 CONTINUE
+ H0 = 0.
+ NRMIN = NRMAX + 1
+ 202 CONTINUE
+* NRMIN = Couronne la Plus Interne
+* NRMIN = Le Plus Petit Rayon a Prendre en Compte
+* DXMIN - H0 = Limite entre Couronnes West et Est
+*
+* Recherche du Premier Rayon
+ H2 = ZHW * ZHW
+ R2 = D2 + H2
+ RW = SQRT(R2)
+*
+* Pente pour le Rayon West
+* Pente pour le Cote
+* PR = - DELTAC / ZHW
+* PC = PENTES(ISXW)
+*
+* -------------
+* Debut
+* Couronne West
+* -------------
+*
+* Aucune Couronne
+ IF (NRMIN .GT. NRMAX) THEN
+ ZHR = ZHE + DXMAX
+ ICC = NRMIN
+*
+* Couronne a Gauche (ZHW <0)
+ ELSEIF (ZHW .LE. (- 0.5 * H0)) THEN
+ IF ((- ZHW * PENTES(ISXW)) .GT. DR) THEN
+ RW = RW - DXMIN
+ ELSE
+ RW = RW + DXMIN
+ ENDIF
+ DO 221 ICC = NRMAX, NRMIN, -1
+ ZRR = RAYONS(ICC)
+ IF (ZRR .LE. RW) THEN
+ GOTO 222
+ ENDIF
+ 221 CONTINUE
+ ICC = NRMIN - 1
+ 222 CONTINUE
+ IF (ICC .LT. NRMIN) THEN
+ ZHR = H0
+ ELSE
+ R2 = ZRR * ZRR
+ H2 = R2 - D2
+ ZHR = - SQRT(H2)
+ ENDIF
+ ICC = ICC + 1
+*
+* Couronne Centrale
+ ELSEIF (ZHW .LE. (0.5 * H0)) THEN
+ ICC = NRMIN
+ ZHR = H0
+*
+* Couronne a Droite (ZHW >0)
+ ELSE
+ IF ((ZHW * PENTES(ISXW)) .GT. - DR) THEN
+ RW = RW + DXMIN
+ ELSE
+ RW = RW - DXMIN
+ ENDIF
+ DO 231 ICC = NRMIN, NRMAX
+ ZRR = RAYONS(ICC)
+ IF (ZRR .GE. RW) THEN
+ GOTO 232
+ ENDIF
+ 231 CONTINUE
+ ICC = NRMAX + 1
+ 232 CONTINUE
+ IF (ICC .GT. NRMAX) THEN
+ ZHR = ZHE + DXMAX
+ ELSE
+ R2 = ZRR * ZRR
+ H2 = R2 - D2
+ ZHR = SQRT(H2)
+ ENDIF
+*
+ ENDIF
+* -------------
+* Couronne West
+* Fin
+* -------------
+*
+* Premiere Position Courante
+ IHC = 1
+ ISC = ISCW
+ ZHC = ZHW
+*
+ IXRAYO(IHC) = ICC
+ HXRAYO(IHC) = ZHC
+*
+* Cote West Intersecte Rayon
+ IF (ZHR .LE. (ZHC + DXMIN)) THEN
+ IF (ZHR .GT. 0.) THEN
+ ICC = ICC + 1
+ ELSE
+ ICC = ICC - 1
+ ENDIF
+ ZHC = MAX (ZHR, ZHC)
+*
+ IHC = IHC + 1
+ IXRAYO(IHC) = ICC
+ HXRAYO(IHC) = ZHC
+*
+ IF (ICC .GT. NRMAX) THEN
+ ZHR = ZHE + DXMAX
+ ELSEIF (ICC .EQ. NRMIN) THEN
+ ZHR = H0
+ ELSEIF (ZHR .LT. 0.) THEN
+ R2 = RAYONS(ICC-1) * RAYONS(ICC-1)
+ H2 = R2 - D2
+ ZHR = - SQRT(H2)
+ ELSE
+ R2 = RAYONS(ICC) * RAYONS(ICC)
+ H2 = R2 - D2
+ ZHR = SQRT(H2)
+ ENDIF
+ ENDIF
+*
+* Suivant
+ 555 CONTINUE
+ ZHN = MIN (ZHE, ZHR, ZHS)
+*
+* +++
+* Debut
+* point Courant (<Est)
+ IF ((ZHN + DXMIN) .LT. ZHE) THEN
+*
+* Increment Rayon
+* Increment Secteur
+ AJRAYO = (ZHN + DXMIN) .GE. ZHR
+ AJSECT = (ZHN + DXMIN) .GE. ZHS
+*
+* soit Intersection Rayon et Secteur
+* soit Intersection SEULEMENT Rayon
+* soit Intersection SEULEMENT Secteur (non Rayon, non Cote)
+*
+ IF (AJRAYO) THEN
+ IF (AJSECT) THEN
+ NBAJ = 2
+* Choix pentes
+* PR = - DELTAC / ZHN
+* PS = ANGLES(ISC)
+ IF (ZHN .GT. 0.) THEN
+ ZZ = ZHN * ANGLES(ISC)
+ IF (ZZ .LT. - DR) THEN
+ AJOU(1) = 'Sect'
+ AJOU(2) = 'Rayo'
+ ELSE
+ AJOU(1) = 'Rayo'
+ AJOU(2) = 'Sect'
+ ENDIF
+ ELSE
+ ZZ = - ZHN * ANGLES(ISC)
+ IF (ZZ .LT. DR) THEN
+ AJOU(1) = 'Sect'
+ AJOU(2) = 'Rayo'
+ ELSE
+ AJOU(1) = 'Rayo'
+ AJOU(2) = 'Sect'
+ ENDIF
+ ENDIF
+* Autre Cas
+* =soit Intersection SEULEMENT Rayon (non Secteur)
+* =soit Intersection SEULEMENT Secteur (non Rayon, non Cote)
+ ELSE
+ NBAJ = 1
+ AJOU(1) = 'Rayo'
+ ENDIF
+ ELSE
+ NBAJ = 1
+ AJOU(1) = 'Sect'
+ ENDIF
+*
+ DO 321 IAJ = 1, NBAJ
+*
+ IHC = IHC + 1
+*
+* Avance Rayon
+ IF (AJOU(IAJ) .EQ. 'Rayo') THEN
+*
+ ZZ = MAX (ZHR, ZHC)
+ HXRAYO(IHC) = ZZ
+ ZHC = ZZ
+*
+ IF (ZHR .GT. 0.) THEN
+ ICC = ICC + 1
+ IF (ICC .GT. NRMAX) THEN
+ ZHR = ZHE + DXMAX
+ ELSE
+ R2 = RAYONS(ICC) * RAYONS(ICC)
+ H2 = R2 - D2
+ ZHR = SQRT(H2)
+ ENDIF
+ ELSE
+ ICC = ICC - 1
+ IF (ICC .EQ. NRMIN) THEN
+ ZHR = H0
+ ELSE
+ R2 = RAYONS(ICC-1) * RAYONS(ICC-1)
+ H2 = R2 - D2
+ ZHR = - SQRT(H2)
+ ENDIF
+ ENDIF
+ IXRAYO(IHC) = ICC
+*
+* Avance Secteur
+ ELSE
+*
+ ZZ = MAX (ZHS, ZHC)
+ HXRAYO(IHC) = ZZ
+ ZHC = ZZ
+*
+ IXRAYO(IHC) = ICC
+ ISC = ISC + 1
+ IF (ISC .GT. ISCMAX) THEN
+ ZHS = ZHE + DXMAX
+ ELSE
+ ZHS = DELTAC * ANGLES(ISC)
+ ENDIF
+ ENDIF
+ 321 CONTINUE
+*
+ GOTO 555
+*
+ ENDIF
+* point Courant (<Est)
+* Fin
+* +++
+*
+* Debut
+* point Final (Est)
+*
+* Intersection Couronne
+ IF ((ZHN + DXMIN) .GE. ZHR) THEN
+* ? Pentes
+* PR = - DELTAC / ZHR
+* PC = PENTES(ISXE)
+ IF (ZHR .GT. 0.) THEN
+ ZZ = ZHR * PENTES(ISXE)
+ AJRAYO = ZZ .GT. - DR
+ ELSE
+ ZZ = - ZHR * PENTES(ISXE)
+ AJRAYO = ZZ .GT. DR
+ ENDIF
+ IF (AJRAYO) THEN
+ IHC = IHC + 1
+ ZZ = MIN(ZHR, ZHE)
+ ZZ = MAX(ZZ , ZHC)
+ HXRAYO(IHC) = ZZ
+ ZHC = ZZ
+ IF (ZHR .GT. 0.) THEN
+ ICC = ICC + 1
+ ELSE
+ ICC = ICC - 1
+ ENDIF
+ IXRAYO(IHC) = ICC
+*
+ ENDIF
+ ENDIF
+*
+ HXRAYO(IHC+1) = MAX(ZHC, ZHE)
+ NHMAX = IHC
+ ISCE = ISC
+*
+ RETURN
+ END