GEANT321/ggeom/gmediv.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:51  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GMEDIV (JVO, IN, XC, IFL)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *   Updates the common /GCVOLU/ and the structure JGPAR          *
  17 C.    *     for contents defined by division.                          *
  18 C.    *                                                                *
  19 C.    *   For IFL nonzero, it also checks if the point XC is inside    *
  20 C.    *     the content. It returns IN = 0, if the point is outside.   *
  21 C.    *     Otherwise, it transforms XC in the local system.           *
  22 C.    *                                                                *
  23 C.    *   For IFL zero, IN is returned 0, if IN > NDIV.                *
  24 C.    *                                                                *
  25 C.    *   Input : JVO, IN, XC, IFL                                     *
  26 C.    *   Output : IN, XC                                              *
  27 C.    *                                                                *
  28 C.    *   Called by : GDRAW, GMEDIA                                    *
  29 C.    *   Authors   : S.Banerjee, R.Brun, F.Bruyant, A.McPherson       *
  30 C.    *                                                                *
  31 C.    ******************************************************************
  32 C.
  33 #include "geant321/gcbank.inc"
  34 #include "geant321/gconsp.inc"
  35 #include "geant321/gcpoly.inc"
  36 #include "geant321/gcvolu.inc"
  37 #if !defined(CERNLIB_SINGLE)
  38       DOUBLE PRECISION DPHIO, TPIDEG, ONE
  39 #endif
  40       DIMENSION  XC(*)
  41       REAL       X0(3)
  42       INTEGER    IDTYP(3,12)
  43       PARAMETER (TPIDEG=360,ONE=1)
  44       SAVE IDTYP
  45 C.
  46       DATA  IDTYP / 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 2, 3, 1,
  47      +              2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 4, 3, 1, 1, 1,
  48      +              2, 3, 1, 2, 3, 1/
  49 C.
  50 C.    ------------------------------------------------------------------
  51 C.
  52       JDIV  = LQ(JVO-1)
  53       ISH   = Q(JVO+2)
  54       IAXIS = Q(JDIV+1)
  55       IVOT  = Q(JDIV+2)
  56       JVOT  = LQ(JVOLUM-IVOT)
  57       IF (NLEVEL.LT.NLDEV(NLEVEL)) THEN
  58          JPAR = 0
  59       ELSE
  60 *       (case with structure JVOLUM locally developed)
  61          JPAR = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL))))
  62          IF (NLEVEL.EQ.NLDEV(NLEVEL)) GO TO 20
  63          DO 10 ILEV = NLDEV(NLEVEL), NLEVEL-1
  64             IF (IQ(JPAR+1).EQ.0) THEN
  65                JPAR = LQ(JPAR-LINDEX(ILEV+1))
  66                IF (JPAR.EQ.0) GO TO 20
  67             ELSE IF (IQ(JPAR-3).GT.1) THEN
  68                JPAR = LQ(JPAR-LINDEX(ILEV+1))
  69             ELSE
  70                JPAR = LQ(JPAR-1)
  71             ENDIF
  72             IF (ILEV.EQ.NLEVEL-1) THEN
  73                NDIV  = IQ(JPAR+1)
  74                ORIG  =  Q(JPAR+2)
  75                SDIV  =  Q(JPAR+3)
  76             ENDIF
  77    10    CONTINUE
  78          GO TO 30
  79       ENDIF
  80 *      (normal case)
  81    20 NDIV  = Q(JDIV+3)
  82       ORIG  = Q(JDIV+4)
  83       SDIV  = Q(JDIV+5)
  84 *
  85    30 IDT = IDTYP(IAXIS,ISH)
  86       IF (IFL.NE.0) THEN
  87          IF (IDT.EQ.1) THEN
  88 *
  89 *         Division along X, Y or Z axis
  90 *
  91             XTT = XC(IAXIS)
  92             IF (ISH.EQ.10) THEN
  93                IF (IAXIS.NE.3) THEN
  94                   XTT = XTT - Q(LQ(JGPAR-NLEVEL)+IAXIS+4) * XC(3)
  95                   IF (IAXIS.EQ.1) THEN
  96                      YT  = XC(2) - Q(LQ(JGPAR-NLEVEL)+6) * XC(3)
  97                      XTT = XTT - Q(LQ(JGPAR-NLEVEL)+4) * YT
  98                   ENDIF
  99                ENDIF
 100             ENDIF
 101             IN = (XTT -ORIG)/SDIV +1
 102          ELSE IF (IDT.EQ.2) THEN
 103 *
 104 *          Division along R axis
 105 *
 106             R = XC(1)**2 + XC(2)**2
 107             IF (ISH.EQ.9) R = R + XC(3)**2
 108             R = SQRT (R)
 109             IF (ISH.EQ.5.OR.ISH.EQ.6.OR.ISH.EQ.9) THEN
 110                IN = (R - ORIG) / SDIV + 1
 111             ELSE IF (ISH.EQ.7.OR.ISH.EQ.8) THEN
 112                IPAR = LQ(JGPAR-NLEVEL)
 113                DR   = 0.5 * (Q(IPAR+4) - Q(IPAR+2)) / Q(IPAR+1)
 114                RMN  = 0.5 * (Q(IPAR+4) + Q(IPAR+2)) + DR * XC(3)
 115                DR   = 0.5 * (Q(IPAR+5) - Q(IPAR+3)) / Q(IPAR+1)
 116                RMX  = 0.5 * (Q(IPAR+5) + Q(IPAR+3)) + DR * XC(3)
 117                STP  = (RMX - RMN) / NDIV
 118                IN   = (R - RMN) / STP + 1
 119             ELSE
 120                IPAR = LQ(JGPAR-NLEVEL)
 121                IF (ISH.EQ.12) THEN
 122                   IPT = IPAR + 1
 123                ELSE
 124                   IPT = IPAR + 2
 125                ENDIF
 126                IF (IZSEC.GT.0) THEN
 127                   IPT = IPT + 3 * IZSEC
 128                ELSE
 129                   NZ  = Q(IPT+2)
 130                   DO 40 IZ = 1, NZ-1
 131                      IF ((XC(3)-Q(IPT+3*IZ))*(XC(3)-Q(IPT+3*IZ+3))
 132      +               .LE.0.) THEN
 133                         IZSEC = IZ
 134                         IPT = IPT + 3 * IZSEC
 135                         GO TO 50
 136                      ENDIF
 137    40             CONTINUE
 138                   IN  = 0
 139                   GO TO 60
 140                ENDIF
 141    50          POR1 = (Q(IPT+3) - XC(3)) / (Q(IPT+3) - Q(IPT))
 142                POR2 = (XC(3) - Q(IPT)) / (Q(IPT+3) - Q(IPT))
 143                RMN  = Q(IPT+1) * POR1 + Q(IPT+4) * POR2
 144                RMX  = Q(IPT+2) * POR1 + Q(IPT+5) * POR2
 145                IF (ISH.EQ.11) THEN
 146                   NPDV = Q(IPAR+3)
 147                   DPH  = Q(IPAR+2) / NPDV
 148                   IF (IPSEC.LE.0) THEN
 149                      IF (XC(1).NE.0..OR.XC(2).NE.0.) THEN
 150                         PHI  = RADDEG * ATAN2 (XC(2), XC(1))
 151                      ELSE
 152                         PHI  = 0.0
 153                      ENDIF
 154                      PH0 = PHI-Q(IPAR+1)
 155                      SG = SIGN(1.0,PH0)
 156                      PH0 = MOD( ABS(PH0), 360.0 )
 157                      IF(SG.LE.0.0) PH0 = 360.0-PH0
 158                      IPSEC= PH0/DPH + 1
 159                   ENDIF
 160                   PH   = DEGRAD * (Q(IPAR+1) + (IPSEC - 0.5) * DPH)
 161                   R    = XC(1) * COS(PH) + XC(2) * SIN(PH)
 162                ENDIF
 163                STP = (RMX - RMN) / NDIV
 164                IN  = (R - RMN) / STP + 1
 165             ENDIF
 166          ELSE IF (IDT.EQ.3) THEN
 167 *
 168 *          Division along Phi axis
 169 *
 170             IF (XC(1).NE.0..OR.XC(2).NE.0.) THEN
 171                PHI = RADDEG * ATAN2 (XC(2), XC(1))
 172             ELSE
 173                PHI = 0.
 174             ENDIF
 175             DPHIO = PHI-ORIG
 176             SG = SIGN(ONE,DPHIO)
 177             DPHIO = MOD( ABS(DPHIO), TPIDEG)
 178             IF(SG.LE.0.0) DPHIO=TPIDEG-DPHIO
 179             IN = DPHIO/SDIV+1
 180          ELSE IF (IDT.EQ.4) THEN
 181 *
 182 *          Division along Theta axis
 183 *
 184             IF (XC(3).NE.0.0) THEN
 185                RXY  = SQRT (XC(1)**2 + XC(2)**2)
 186                THET = RADDEG * ATAN (RXY/XC(3))
 187                IF (THET.LT.0.0)  THET = THET + 180.0
 188             ELSE
 189                THET = 90.0
 190             ENDIF
 191             IN   = (THET - ORIG) / SDIV + 1
 192          ENDIF
 193       ENDIF
 194 *
 195    60 IF (IN.GT.NDIV) IN = 0
 196       IF (IN.LE.0) GO TO 999
 197 *
 198       IF (JPAR.NE.0) THEN
 199          IF (IQ(JPAR-3).GT.1) THEN
 200             JPAR = LQ(JPAR-IN)
 201          ELSE
 202             JPAR = LQ(JPAR-1)
 203          ENDIF
 204          JPAR = JPAR + 5
 205          NPAR = IQ(JPAR)
 206       ELSE
 207          NPAR = Q(JVOT+5)
 208          JPAR = JVOT + 6
 209       ENDIF
 210 *
 211 *      Volume found at deeper level
 212 *
 213       NL1    = NLEVEL
 214       NLEVEL = NLEVEL +1
 215       LVOLUM(NLEVEL) = IVOT
 216       NAMES(NLEVEL)  = IQ(JVOLUM+IVOT)
 217       NUMBER(NLEVEL) = IN
 218       LINDEX(NLEVEL) = IN
 219       LINMX(NLEVEL)  = NDIV
 220       GONLY(NLEVEL)  = GONLY(NL1)
 221       IF (LQ(LQ(JVOLUM-IVOT)).EQ.0) THEN
 222          NLDEV(NLEVEL) = NLDEV(NL1)
 223       ELSE
 224          NLDEV(NLEVEL) = NLEVEL
 225       ENDIF
 226 *
 227       IF (IDT.EQ.1) THEN
 228          X0(1) = 0.0
 229          X0(2) = 0.0
 230          X0(3) = 0.0
 231          X0(IAXIS) = ORIG + (IN - 0.5) * SDIV
 232          IF (ISH.EQ.4.OR.(ISH.EQ.10.AND.IAXIS.NE.1)) THEN
 233             CALL GCENT (IAXIS, X0)
 234          ENDIF
 235          IF (IFL.NE.0) THEN
 236             XC(1) = XC(1) - X0(1)
 237             XC(2) = XC(2) - X0(2)
 238             XC(3) = XC(3) - X0(3)
 239          ENDIF
 240 C*****  Code Expanded From Routine:  GTRMUL
 241 C.
 242 C.    ------------------------------------------------------------------
 243 C.
 244          IF (GRMAT(10,NL1) .EQ. 0.0) THEN
 245             GTRAN(1,NLEVEL) = GTRAN(1,NL1) + X0(1)
 246             GTRAN(2,NLEVEL) = GTRAN(2,NL1) + X0(2)
 247             GTRAN(3,NLEVEL) = GTRAN(3,NL1) + X0(3)
 248             DO 70 I = 1, 10, 2
 249                GRMAT(I,NLEVEL) = GRMAT(I,NL1)
 250                GRMAT(I+1,NLEVEL) = GRMAT(I+1,NL1)
 251    70       CONTINUE
 252          ELSE
 253 C
 254             DXTEM1 = X0(1)*GRMAT(1,NL1) + X0(2)*GRMAT(4,NL1) + X0(3)*
 255      +      GRMAT( 7,NL1)
 256             DXTEM2 = X0(1)*GRMAT(2,NL1) + X0(2)*GRMAT(5,NL1) + X0(3)*
 257      +      GRMAT( 8,NL1)
 258             DXTEM3 = X0(1)*GRMAT(3,NL1) + X0(2)*GRMAT(6,NL1) + X0(3)*
 259      +      GRMAT( 9,NL1)
 260             DO 80 I = 1, 10, 2
 261                GRMAT(I,NLEVEL) = GRMAT(I,NL1)
 262                GRMAT(I+1,NLEVEL) = GRMAT(I+1,NL1)
 263    80       CONTINUE
 264             GTRAN(1,NLEVEL) = GTRAN(1,NL1) + DXTEM1
 265             GTRAN(2,NLEVEL) = GTRAN(2,NL1) + DXTEM2
 266             GTRAN(3,NLEVEL) = GTRAN(3,NL1) + DXTEM3
 267          ENDIF
 268 C*****  End of Code Expanded From Routine:  GTRMUL
 269 *
 270       ELSE IF (IDT.EQ.3.OR.IDT.EQ.4) THEN
 271          IF (IDT.EQ.3) THEN
 272             PH0  = DEGRAD * (ORIG + (IN - 0.5) * SDIV)
 273             CPHR = COS (PH0)
 274             SPHR = SIN (PH0)
 275          ELSE
 276             PH0  = 0.0
 277             CPHR = 1.0
 278             SPHR = 0.0
 279          ENDIF
 280          GTRAN(1,NLEVEL) = GTRAN(1,NL1)
 281          GRMAT(1,NLEVEL) = GRMAT(1,NL1)*CPHR + GRMAT(4,NL1)*SPHR
 282          GRMAT(4,NLEVEL) = GRMAT(4,NL1)*CPHR - GRMAT(1,NL1)*SPHR
 283          GRMAT(7,NLEVEL) = GRMAT(7,NL1)
 284          GTRAN(2,NLEVEL) = GTRAN(2,NL1)
 285          GRMAT(2,NLEVEL) = GRMAT(2,NL1)*CPHR + GRMAT(5,NL1)*SPHR
 286          GRMAT(5,NLEVEL) = GRMAT(5,NL1)*CPHR - GRMAT(2,NL1)*SPHR
 287          GRMAT(8,NLEVEL) = GRMAT(8,NL1)
 288          GTRAN(3,NLEVEL) = GTRAN(3,NL1)
 289          GRMAT(3,NLEVEL) = GRMAT(3,NL1)*CPHR + GRMAT(6,NL1)*SPHR
 290          GRMAT(6,NLEVEL) = GRMAT(6,NL1)*CPHR - GRMAT(3,NL1)*SPHR
 291          GRMAT(9,NLEVEL) = GRMAT(9,NL1)
 292          IF (IFL.NE.0) THEN
 293             XTT   = XC(1) * CPHR + XC(2) * SPHR
 294             XC(2) = XC(2) * CPHR - XC(1) * SPHR
 295             XC(1) = XTT
 296          ENDIF
 297          IF (PH0.EQ.0.0.AND.GRMAT(10,NL1).EQ.0.0) THEN
 298             GRMAT(10,NLEVEL) = 0.0
 299          ELSE
 300             GRMAT(10,NLEVEL) = 1.0
 301          ENDIF
 302          IF (ISH.EQ.11) IPSEC = 1
 303 *
 304       ELSE
 305          GTRAN(1,NLEVEL) = GTRAN(1,NL1)
 306          GTRAN(2,NLEVEL) = GTRAN(2,NL1)
 307          GTRAN(3,NLEVEL) = GTRAN(3,NL1)
 308          DO 90 I = 1, 10, 2
 309             GRMAT(I,NLEVEL) = GRMAT(I,NL1)
 310             GRMAT(I+1,NLEVEL) = GRMAT(I+1,NL1)
 311    90    CONTINUE
 312       ENDIF
 313 *
 314       IQ(JGPAR+NLEVEL) = NPAR
 315       LQ(JGPAR-NLEVEL) = JPAR
 316 *                                                             END GMEDIV
 317   999 END