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