]> git.uio.no Git - u/mrichter/AliRoot.git/blob - GEANT321/peanut/sbcomp.F
Initial version
[u/mrichter/AliRoot.git] / GEANT321 / peanut / sbcomp.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:22:03  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.46  by  S.Giani
11 *-- Author :
12 *$ CREATE SBCOMP.FOR
13 *COPY SBCOMP
14 *
15 *=== sbcomp ===========================================================*
16 *
17       SUBROUTINE SBCOMP ( SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1,
18      &                    SBHAL1 )
19  
20 #include "geant321/dblprc.inc"
21 #include "geant321/dimpar.inc"
22 #include "geant321/iounit.inc"
23 *
24 *----------------------------------------------------------------------*
25 *----------------------------------------------------------------------*
26 *
27 #include "geant321/nucgeo.inc"
28 *
29       SAVE  R0HLA , R0HLB , R0SKA , R0SKB , R0CEA , R0CEB ,
30      &      R1CEA , R1CEB , R1SKA , R1SKB , R1HLA , R1HLB ,
31      &      SQRH0A, SQRH0B, SQRS0A, SQRS0B, SQRC0A, SQRC0B,
32      &      SQRC1A, SQRC1B, SQRS1A, SQRS1B, SQRH1A, SQRH1B,
33      &      BIMPSQ
34 #include "geant321/nucstf.inc"
35 *
36       IF ( ABS (R0TRAJ) .LT. BIMPCT .OR. ABS (R1TRAJ) .LT. BIMPCT )
37      &   STOP 'BIMPCT'
38       BIMPSQ = BIMPCT * BIMPCT
39       IF ( R0TRAJ .GE. 0.D+00 ) THEN
40          SBHAL0 = 0.D+00
41          SBSKI0 = 0.D+00
42          SBCEN0 = 0.D+00
43       ELSE
44          IF ( BIMPCT .GE. RADIU1 ) THEN
45             SBSKI0 = 0.D+00
46             SBCEN0 = 0.D+00
47             SBCEN1 = 0.D+00
48             SBSKI1 = 0.D+00
49             R0HLA = ABS (R0TRAJ)
50             R0HLB = MAX ( BIMPCT, - R1TRAJ )
51             SQRH0A = SQRT ( ( R0HLA - BIMPCT ) * ( R0HLA + BIMPCT ) )
52             SQRH0B = SQRT ( ( R0HLB - BIMPCT ) * ( R0HLB + BIMPCT ) )
53             SBHAL0 = RHOSKN * ( ( 1.D+00 + RADIU1 / HALODP ) * ( SQRH0A
54      &             - SQRH0B ) - 0.5D+00 * ( R0HLA / HALODP * SQRH0A
55      &             - R0HLB / HALODP * SQRH0B + BIMPSQ / HALODP
56      &             * LOG ( ( R0HLA + SQRH0A ) / ( R0HLB + SQRH0B ) ) ) )
57          ELSE IF ( R0TRAJ .GE. -RADIU0 ) THEN
58             SBHAL0 = 0.D+00
59             SBSKI0 = 0.D+00
60             R0CEA = - R0TRAJ
61             R0CEB = MAX ( BIMPCT, - R1TRAJ )
62             SQRC0A = SQRT ( ( R0CEA - BIMPCT ) * ( R0CEA + BIMPCT ) )
63             SQRC0B = SQRT ( ( R0CEB - BIMPCT ) * ( R0CEB + BIMPCT ) )
64             SBCEN0 = RHOCEN * ( SQRC0A - SQRC0B )
65          ELSE IF ( R0TRAJ .GE. -RADIU1 ) THEN
66             SBHAL0 = 0.D+00
67             RHELP = MAX ( BIMPCT, - R1TRAJ )
68             R0SKA = ABS (R0TRAJ)
69             R0SKB = MAX ( RHELP, RADIU0 )
70             SQRS0A = SQRT ( ( R0SKA - BIMPCT ) * ( R0SKA + BIMPCT ) )
71             SQRS0B = SQRT ( ( R0SKB - BIMPCT ) * ( R0SKB + BIMPCT ) )
72             SBSKI0 = RHOCEN * ( ( 1.D+00 + RADIU0 / SKNEFF ) * ( SQRS0A
73      &             - SQRS0B ) - 0.5D+00 * ( R0SKA / SKNEFF * SQRS0A
74      &             - R0SKB / SKNEFF * SQRS0B + BIMPSQ / SKNEFF
75      &             * LOG ( ( R0SKA + SQRS0A ) / ( R0SKB + SQRS0B ) ) ) )
76             IF ( RHELP .LT. RADIU0 ) THEN
77                R0CEA = RADIU0
78                R0CEB = RHELP
79                SQRC0A = SQRT ( ( R0CEA - BIMPCT ) * ( R0CEA + BIMPCT ) )
80                SQRC0B = SQRT ( ( R0CEB - BIMPCT ) * ( R0CEB + BIMPCT ) )
81                SBCEN0 = RHOCEN * ( SQRC0A - SQRC0B )
82             ELSE
83                SBCEN0 = 0.D+00
84             END IF
85          ELSE
86             RHELP = MAX ( BIMPCT, - R1TRAJ )
87             R0HLA = ABS (R0TRAJ)
88             R0HLB = MAX ( RHELP, RADIU1 )
89             SQRH0A = SQRT ( ( R0HLA - BIMPCT ) * ( R0HLA + BIMPCT ) )
90             SQRH0B = SQRT ( ( R0HLB - BIMPCT ) * ( R0HLB + BIMPCT ) )
91             SBHAL0 = RHOSKN * ( ( 1.D+00 + RADIU1 / HALODP ) * ( SQRH0A
92      &             - SQRH0B ) - 0.5D+00 * ( R0HLA / HALODP * SQRH0A
93      &             - R0HLB / HALODP * SQRH0B + BIMPSQ / HALODP
94      &             * LOG ( ( R0HLA + SQRH0A ) / ( R0HLB + SQRH0B ) ) ) )
95             IF ( RHELP .LT. RADIU1 ) THEN
96                R0SKA = RADIU1
97                R0SKB = MAX ( RHELP, RADIU0 )
98                SQRS0A = SQRT ( ( R0SKA - BIMPCT ) * ( R0SKA + BIMPCT ) )
99                SQRS0B = SQRT ( ( R0SKB - BIMPCT ) * ( R0SKB + BIMPCT ) )
100                SBSKI0 = RHOCEN * ( ( 1.D+00 + RADIU0 / SKNEFF )
101      &                * ( SQRS0A - SQRS0B ) - 0.5D+00 * ( R0SKA / SKNEFF
102      &                * SQRS0A - R0SKB / SKNEFF * SQRS0B + BIMPSQ
103      &                / SKNEFF * LOG ( ( R0SKA + SQRS0A ) / ( R0SKB
104      &                + SQRS0B ) ) ) )
105                IF ( RHELP .LT. RADIU0 ) THEN
106                   R0CEA = RADIU0
107                   R0CEB = MAX ( BIMPCT, - R1TRAJ )
108                   SQRC0A = SQRS0B
109                   SQRC0B = SQRT ( ( R0CEB - BIMPCT )*( R0CEB + BIMPCT ))
110                   SBCEN0 = RHOCEN * ( SQRC0A - SQRC0B )
111                ELSE
112                   SBCEN0 = 0.D+00
113                END IF
114             ELSE
115                SBSKI0 = 0.D+00
116                SBCEN0 = 0.D+00
117             END IF
118          END IF
119       END IF
120       IF ( R1TRAJ .EQ. - R0TRAJ ) THEN
121          R1HLA  = R0HLB
122          R1HLB  = R0HLA
123          SQRH1A = SQRH0B
124          SQRH1B = SQRH0A
125          R1SKA  = R0SKB
126          R1SKB  = R0SKA
127          SQRS1A = SQRS0B
128          SQRS1B = SQRS0A
129          R1CEA  = R0CEB
130          R1CEB  = R0CEA
131          SQRC1A = SQRC0B
132          SQRC1B = SQRC0A
133          SBCEN1 = SBCEN0
134          SBSKI1 = SBSKI0
135          SBHAL1 = SBHAL0
136       ELSE IF ( R1TRAJ .LE. 0.D+00 ) THEN
137          SBCEN1 = 0.D+00
138          SBSKI1 = 0.D+00
139          SBHAL1 = 0.D+00
140       ELSE
141          IF ( BIMPCT .GE. RADIU1 ) THEN
142             SBSKI1 = 0.D+00
143             SBCEN1 = 0.D+00
144             R1HLB = R1TRAJ
145             R1HLA = MAX ( BIMPCT, R0TRAJ )
146             SQRH1A = SQRT ( ( R1HLA - BIMPCT ) * ( R1HLA + BIMPCT ) )
147             SQRH1B = SQRT ( ( R1HLB - BIMPCT ) * ( R1HLB + BIMPCT ) )
148             SBHAL1 = RHOSKN * ( ( 1.D+00 + RADIU1 / HALODP ) * ( SQRH1B
149      &             - SQRH1A ) - 0.5D+00 * ( R1HLB / HALODP * SQRH1B
150      &             - R1HLA / HALODP * SQRH1A + BIMPSQ / HALODP
151      &             * LOG ( ( R1HLB + SQRH1B ) / ( R1HLA + SQRH1A ) ) ) )
152          ELSE IF ( R1TRAJ .LE. RADIU0 ) THEN
153             SBHAL1 = 0.D+00
154             SBSKI1 = 0.D+00
155             R1CEB = R1TRAJ
156             R1CEA = MAX ( BIMPCT, R0TRAJ )
157             SQRC1A = SQRT ( ( R1CEA - BIMPCT ) * ( R1CEA + BIMPCT ) )
158             SQRC1B = SQRT ( ( R1CEB - BIMPCT ) * ( R1CEB + BIMPCT ) )
159             SBCEN1 = RHOCEN * ( SQRC1B - SQRC1A )
160          ELSE IF ( R1TRAJ .LE. RADIU1 ) THEN
161             SBHAL1 = 0.D+00
162             RHELP = MAX ( BIMPCT, R0TRAJ )
163             R1SKB = R1TRAJ
164             R1SKA = MAX ( RHELP, RADIU0 )
165             SQRS1A = SQRT ( ( R1SKA - BIMPCT ) * ( R1SKA + BIMPCT ) )
166             SQRS1B = SQRT ( ( R1SKB - BIMPCT ) * ( R1SKB + BIMPCT ) )
167             SBSKI1 = RHOCEN * ( ( 1.D+00 + RADIU0 / SKNEFF ) * ( SQRS1B
168      &             - SQRS1A ) - 0.5D+00 * ( R1SKB / SKNEFF * SQRS1B
169      &             - R1SKA / SKNEFF * SQRS1A + BIMPSQ / SKNEFF
170      &             * LOG ( ( R1SKB + SQRS1B ) / ( R1SKA + SQRS1A ) ) ) )
171             IF ( RHELP .LT. RADIU0 ) THEN
172                R1CEB = RADIU0
173                R1CEA = RHELP
174                SQRC1A = SQRT ( ( R1CEA - BIMPCT ) * ( R1CEA + BIMPCT ) )
175                SQRC1B = SQRT ( ( R1CEB - BIMPCT ) * ( R1CEB + BIMPCT ) )
176                SBCEN1 = RHOCEN * ( SQRC1B - SQRC1A )
177             ELSE
178                SBCEN1 = 0.D+00
179             END IF
180          ELSE
181             RHELP = MAX ( BIMPCT, R0TRAJ )
182             R1HLB = R1TRAJ
183             R1HLA = MAX ( RHELP, RADIU1 )
184             SQRH1A = SQRT ( ( R1HLA - BIMPCT ) * ( R1HLA + BIMPCT ) )
185             SQRH1B = SQRT ( ( R1HLB - BIMPCT ) * ( R1HLB + BIMPCT ) )
186             SBHAL1 = RHOSKN * ( ( 1.D+00 + RADIU1 / HALODP ) * ( SQRH1B
187      &             - SQRH1A ) - 0.5D+00 * ( R1HLB / HALODP * SQRH1B
188      &             - R1HLA / HALODP * SQRH1A + BIMPSQ / HALODP
189      &             * LOG ( ( R1HLB + SQRH1B ) / ( R1HLA + SQRH1A ) ) ) )
190             IF ( RHELP .LT. RADIU1 ) THEN
191                R1SKB = RADIU1
192                R1SKA = MAX ( RHELP, RADIU0 )
193                SQRS1A = SQRT ( ( R1SKA - BIMPCT ) * ( R1SKA + BIMPCT ) )
194                SQRS1B = SQRT ( ( R1SKB - BIMPCT ) * ( R1SKB + BIMPCT ) )
195                SBSKI1 = RHOCEN * ( ( 1.D+00 + RADIU0 / SKNEFF )
196      &                * ( SQRS1B - SQRS1A ) - 0.5D+00 * ( R1SKB / SKNEFF
197      &                * SQRS1B - R1SKA / SKNEFF * SQRS1A + BIMPSQ
198      &                / SKNEFF * LOG ( ( R1SKB + SQRS1B ) / ( R1SKA
199      &                + SQRS1A ) ) ) )
200                IF ( RHELP .LT. RADIU0 ) THEN
201                   R1CEB = RADIU0
202                   R1CEA = RHELP
203                   SQRC1B = SQRS1A
204                   SQRC1A = SQRT ( ( R1CEA - BIMPCT )*( R1CEA + BIMPCT ))
205                   SBCEN1 = RHOCEN * ( SQRC1B - SQRC1A )
206                ELSE
207                   SBCEN1 = 0.D+00
208                END IF
209             ELSE
210                SBCEN1 = 0.D+00
211                SBSKI1 = 0.D+00
212             END IF
213          END IF
214       END IF
215       RETURN
216 *
217 *----------------------------------------------------------------------*
218 *----------------------------------------------------------------------*
219 *
220       ENTRY RSCOMP ( SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1 )
221       SBTTOT = SBHAL0 + SBSKI0 + SBCEN0 + SBCEN1 + SBSKI1 + SBHAL1
222       IF ( SBUSED .GT. SBTTOT ) STOP 'RSCOMP'
223       IF ( SBUSED .LT. SBHAL0 ) THEN
224          SBUSEF = SBUSED
225          RDLESS = R0HLA
226          RDMORE = R0HLB
227          SBLESS = 0.D+00
228          SBMORE = SBHAL0
229          DO 1000 IBIS = 1,4
230             RIMPCT = 0.5D+00 * ( RDMORE + RDLESS )
231             SQRIMP = SQRT ( ( RIMPCT - BIMPCT ) * ( RIMPCT + BIMPCT ) )
232             SBRIMP = RHOSKN * ( ( 1.D+00 + RADIU1 / HALODP ) * ( SQRH0A
233      &             - SQRIMP ) - 0.5D+00 * ( R0HLA / HALODP * SQRH0A
234      &             - RIMPCT / HALODP * SQRIMP + BIMPSQ / HALODP
235      &             * LOG ( ( R0HLA + SQRH0A ) / ( RIMPCT + SQRIMP ) ) ))
236             IF ( SBRIMP .GE. SBUSEF ) THEN
237                SBMORE = SBRIMP
238                RDMORE = RIMPCT
239             ELSE
240                SBLESS = SBRIMP
241                RDLESS = RIMPCT
242             END IF
243  1000    CONTINUE
244          RIMPCT = RDLESS + ( SBUSEF - SBLESS ) / ( SBMORE - SBLESS )
245      &          * ( RDMORE - RDLESS )
246          RIMPCT = - RIMPCT
247          SQRIMP = SQRT ( ( RIMPCT - BIMPCT ) * ( RIMPCT + BIMPCT ) )
248       ELSE IF ( SBUSED .LT. SBHAL0 + SBSKI0 ) THEN
249          SBUSEF = SBUSED - SBHAL0
250          RDLESS = R0SKA
251          RDMORE = R0SKB
252          SBLESS = 0.D+00
253          SBMORE = SBSKI0
254          DO 2000 IBIS = 1,4
255             RIMPCT = 0.5D+00 * ( RDMORE + RDLESS )
256             SQRIMP = SQRT ( ( RIMPCT - BIMPCT ) * ( RIMPCT + BIMPCT ) )
257             SBRIMP = RHOCEN * ( ( 1.D+00 + RADIU0 / SKNEFF ) * ( SQRS0A
258      &             - SQRIMP ) - 0.5D+00 * ( R0SKA / SKNEFF * SQRS0A
259      &             - RIMPCT / SKNEFF * SQRIMP + BIMPSQ / SKNEFF
260      &             * LOG ( ( R0SKA + SQRS0A ) / ( RIMPCT + SQRIMP ) ) ))
261             IF ( SBRIMP .GE. SBUSEF ) THEN
262                SBMORE = SBRIMP
263                RDMORE = RIMPCT
264             ELSE
265                SBLESS = SBRIMP
266                RDLESS = RIMPCT
267             END IF
268  2000    CONTINUE
269          RIMPCT = RDLESS + ( SBUSEF - SBLESS ) / ( SBMORE - SBLESS )
270      &          * ( RDMORE - RDLESS )
271          RIMPCT = - RIMPCT
272          SQRIMP = SQRT ( ( RIMPCT - BIMPCT ) * ( RIMPCT + BIMPCT ) )
273       ELSE IF ( SBUSED .LT. SBHAL0 + SBSKI0 + SBCEN0 ) THEN
274          SBUSEF = SBUSED - SBHAL0 - SBSKI0
275          SQRIMP = SQRC0A - SBUSEF / RHOCEN
276          RIMPCT = SQRT ( SQRIMP**2 + BIMPSQ )
277          RIMPCT = - RIMPCT
278       ELSE IF ( SBUSED .LT. SBTTOT - SBSKI1 - SBHAL1 ) THEN
279          SBUSEF = SBUSED - SBHAL0 - SBSKI0 - SBCEN0
280          SQRIMP = SBUSEF / RHOCEN + SQRC1A
281          RIMPCT = SQRT ( SQRIMP**2 + BIMPSQ )
282       ELSE IF ( SBUSED .LT. SBTTOT - SBHAL1 ) THEN
283          SBUSEF = SBUSED - SBHAL0 - SBSKI0 - SBCEN0 - SBCEN1
284          RDLESS = R1SKA
285          RDMORE = R1SKB
286          SBLESS = 0.D+00
287          SBMORE = SBSKI1
288          DO 3000 IBIS = 1,4
289             RIMPCT = 0.5D+00 * ( RDMORE + RDLESS )
290             SQRIMP = SQRT ( ( RIMPCT - BIMPCT ) * ( RIMPCT + BIMPCT ) )
291             SBRIMP = RHOCEN * ( ( 1.D+00 + RADIU0 / SKNEFF ) * ( SQRIMP
292      &             - SQRS1A ) - 0.5D+00 * ( RIMPCT / SKNEFF * SQRIMP
293      &             - R1SKA / SKNEFF * SQRS1A + BIMPSQ / SKNEFF
294      &             * LOG ( ( RIMPCT + SQRIMP ) / ( R1SKA + SQRS1A ) ) ))
295             IF ( SBRIMP .GE. SBUSEF ) THEN
296                SBMORE = SBRIMP
297                RDMORE = RIMPCT
298             ELSE
299                SBLESS = SBRIMP
300                RDLESS = RIMPCT
301             END IF
302  3000    CONTINUE
303          RIMPCT = RDLESS + ( SBUSEF - SBLESS ) / ( SBMORE - SBLESS )
304      &          * ( RDMORE - RDLESS )
305          SQRIMP = SQRT ( ( RIMPCT - BIMPCT ) * ( RIMPCT + BIMPCT ) )
306       ELSE
307          SBUSEF = SBUSED - SBHAL0 - SBSKI0 - SBCEN0 - SBCEN1 - SBSKI1
308          RDLESS = R1HLA
309          RDMORE = R1HLB
310          SBLESS = 0.D+00
311          SBMORE = SBHAL1
312          DO 4000 IBIS = 1,4
313             RIMPCT = 0.5D+00 * ( RDMORE + RDLESS )
314             SQRIMP = SQRT ( ( RIMPCT - BIMPCT ) * ( RIMPCT + BIMPCT ) )
315             SBRIMP = RHOSKN * ( ( 1.D+00 + RADIU1 / HALODP ) * ( SQRIMP
316      &             - SQRH1A ) - 0.5D+00 * ( RIMPCT / HALODP * SQRIMP
317      &             - R1HLA / HALODP * SQRH1A + BIMPSQ / HALODP
318      &             * LOG ( ( RIMPCT + SQRIMP ) / ( R1HLA + SQRH1A ) ) ))
319             IF ( SBRIMP .GE. SBUSEF ) THEN
320                SBMORE = SBRIMP
321                RDMORE = RIMPCT
322             ELSE
323                SBLESS = SBRIMP
324                RDLESS = RIMPCT
325             END IF
326  4000    CONTINUE
327          RIMPCT = RDLESS + ( SBUSEF - SBLESS ) / ( SBMORE - SBLESS )
328      &          * ( RDMORE - RDLESS )
329          SQRIMP = SQRT ( ( RIMPCT - BIMPCT ) * ( RIMPCT + BIMPCT ) )
330       END IF
331       IF ( RIMPCT .GT. 0.D+00 ) THEN
332          XIMPCT = XBIMPC + CXIMPC * SQRIMP
333          YIMPCT = YBIMPC + CYIMPC * SQRIMP
334          ZIMPCT = ZBIMPC + CZIMPC * SQRIMP
335       ELSE
336          XIMPCT = XBIMPC - CXIMPC * SQRIMP
337          YIMPCT = YBIMPC - CYIMPC * SQRIMP
338          ZIMPCT = ZBIMPC - CZIMPC * SQRIMP
339       END IF
340       IF ( BIMPCT .GT. ANGLGB ) THEN
341          TRUFAC = BIMPTR / BIMPCT
342          XIMPTR = XIMPCT * TRUFAC
343          YIMPTR = YIMPCT * TRUFAC
344          ZIMPTR = ZIMPCT * TRUFAC
345          RIMPTR = RIMPCT * TRUFAC
346       ELSE
347          XIMPTR = XIMPCT
348          YIMPTR = YIMPCT
349          ZIMPTR = ZIMPCT
350          RIMPTR = RIMPCT
351       END IF
352       RETURN
353 *=== End of subroutine Sbcomp =========================================*
354       END