summaryrefslogtreecommitdiff
path: root/Dragon/src/SYB7TC.f
blob: a4eca3fbce2df573b6200102bf8b985405f7d3e0 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
*DECK SYB7TC
      SUBROUTINE SYB7TC (DELTAR, DDELTA, ANGLES, NHMAX, IXRAYO,
     &                   NRI,    RAYONS, ZZW, ZZE, ZZR, HXRAYO)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the intersection lenghts of a track in a sectorized hexagonal
* 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
* DELTAR  used to compute the intersection.
* DDELTA  used to compute the intersection.
* ANGLES  angular values (begin at 1 or 2):
*         1=   $\\tan$(-$\\pi$/6-PHI);
*         2=   $\\tan$(+$\\pi$/6-PHI);
*         3=   $\\tan$(3$\\pi$/6-PHI).
* NHMAX   number of intersections.
* IXRAYO  tube indices.
* NRI     number of radii (= NRD-1)
* RAYONS  radius of each cylinder.
* ZZW     position of the west intersection (left).
* ZZE     position of the east intersection (right).
*
*Parameters: output
* ZZR     intersection lenghts.
*
*Parameters: input/output
* HXRAYO  preceding/next intersection lenghts.
*
*-----------------------------------------------------------------------
*
*----
*  SUBROUTINE ARGUMENTS
*----
      INTEGER NHMAX,IXRAYO(NHMAX),NRI
      REAL    DELTAR,DDELTA,ANGLES(3),RAYONS(NRI),ZZW,ZZE,ZZR(NHMAX),
     &        HXRAYO(NHMAX+1)
*
      DDELT2 = DDELTA * DDELTA
      DELTA2 = DELTAR * DELTAR
      GC  = 0.
      HC  = ZZW
      IRC = IXRAYO(1)
      ISC = 0
      HPRAYO = HXRAYO(1)
      HXRAYO(1) = HC
*
      DO IH = 1, NHMAX-1
        GP  = GC
        GC  = 0.
        HP  = HC
        IRP = IRC
        IRR = 0
*
        IRC = IXRAYO(IH+1)
        IF (IRP .EQ. IRC) THEN
          ISC = ISC + 1
          HC = DELTAR * ANGLES(ISC)
        ELSE
          IRR = MIN(IRP, IRC)
*
* Distance
          H2 = RAYONS(IRR) * RAYONS(IRR) - DELTA2
          IF (H2 .GT. 0) THEN
            HC = SQRT(H2)
            IF (IRC .EQ. IRR) THEN
              HC = - HC
            ENDIF
          ELSE
            HC = 0.
          ENDIF
*
        ENDIF
*
* Protection contre les longueurs negatives
        IF (HC .LT. HP) THEN
          HC = HP
          ZZH = 0.
        ELSE
          ZZH = HC - HP
        ENDIF
        ZZH = (ZZH + HXRAYO(IH+1) - HPRAYO) * 0.5
        HPRAYO = HXRAYO(IH+1)
        HXRAYO(IH+1) = HC
*
* Ajout de la courbure
        IF (IRP .NE. IRC) THEN
          H2 = HC - HPRAYO
          H2 = H2 * H2
          XCORDE = H2 + DDELT2
          IF (XCORDE .GT. 0.) THEN
* Surface entre la corde et l'arc
            XUNITE = SQRT(XCORDE) / RAYONS(IRR) / 2.
            XALPHA = ASIN(XUNITE)
            XUNITE = XALPHA - COS(XALPHA) * XUNITE
            GC = XUNITE * RAYONS(IRR) * RAYONS(IRR) / DDELTA
          ELSE
            GC = 0.
          ENDIF
*
          IF (IRC .EQ. IRR) THEN
            GC = - GC
          ENDIF
        ENDIF
*
* Longueur Moyenne
        ZZR(IH) = ZZH + GC - GP
*
      ENDDO
*
* Dernier
        IF (ZZE .LT. HC) THEN
          ZZH = 0.
        ELSE
          ZZH = ZZE - HC
        ENDIF
        ZZH = (ZZH + HXRAYO(NHMAX+1) - HPRAYO) * 0.5
        ZZR(NHMAX) = ZZH - GC
        HXRAYO(NHMAX+1) = ZZE
*
      RETURN
      END