GEANT321/gtrak/grkuta.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:42  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.23  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GRKUTA (CHARGE,STEP,VECT,VOUT)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *  Runge-Kutta method for tracking a particle through a magnetic *
  17 C.    *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
  18 C.    *  Standards, procedure 25.5.20)                                 *
  19 C.    *                                                                *
  20 C.    *  Input parameters                                              *
  21 C.    *       CHARGE    Particle charge                                *
  22 C.    *       STEP      Step size                                      *
  23 C.    *       VECT      Initial co-ords,direction cosines,momentum     *
  24 C.    *  Output parameters                                             *
  25 C.    *       VOUT      Output co-ords,direction cosines,momentum      *
  26 C.    *  User routine called                                           *
  27 C.    *       CALL GUFLD(X,F)                                          *
  28 C.    *                                                                *
  29 C.    *    ==>Called by : <USER>, GUSWIM                               *
  30 C.    *       Authors    R.Brun, M.Hansroul  *********                 *
  31 C.    *                  V.Perevoztchikov (CUT STEP implementation)    *
  32 C.    *                                                                *
  33 C.    *                                                                *
  34 C.    ******************************************************************
  35 C.
  36 #if !defined(CERNLIB_SINGLE)
  37       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  38 #endif
  39       REAL CHARGE, STEP, VECT(*), VOUT(*), F(4)
  40       REAL XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
  41       DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
  42       EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
  43      +            (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
  44 *
  45       PARAMETER (MAXIT = 1992, MAXCUT = 11)
  46       PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
  47       PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
  48       PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
  49       PARAMETER (PISQUA=.986960440109D+01)
  50       PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
  51 *.
  52 *.    ------------------------------------------------------------------
  53 *.
  54 *             This constant is for units CM,GEV/C and KGAUSS
  55 *
  56       ITER = 0
  57       NCUT = 0
  58       DO 10 J=1,7
  59          VOUT(J)=VECT(J)
  60    10 CONTINUE
  61       PINV   = EC * CHARGE / VECT(7)
  62       TL = 0.
  63       H      = STEP
  64 *
  65 *
  66    20 REST  = STEP-TL
  67       IF (ABS(H).GT.ABS(REST)) H = REST
  68       CALL GUFLD(VOUT,F)
  69 *
  70 *             Start of integration
  71 *
  72       X      = VOUT(1)
  73       Y      = VOUT(2)
  74       Z      = VOUT(3)
  75       A      = VOUT(4)
  76       B      = VOUT(5)
  77       C      = VOUT(6)
  78 *
  79       H2     = HALF * H
  80       H4     = HALF * H2
  81       PH     = PINV * H
  82       PH2    = HALF * PH
  83       SECXS(1) = (B * F(3) - C * F(2)) * PH2
  84       SECYS(1) = (C * F(1) - A * F(3)) * PH2
  85       SECZS(1) = (A * F(2) - B * F(1)) * PH2
  86       ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
  87       IF (ANG2.GT.PISQUA) GO TO 40
  88       DXT    = H2 * A + H4 * SECXS(1)
  89       DYT    = H2 * B + H4 * SECYS(1)
  90       DZT    = H2 * C + H4 * SECZS(1)
  91       XT     = X + DXT
  92       YT     = Y + DYT
  93       ZT     = Z + DZT
  94 *
  95 *              Second intermediate point
  96 *
  97       EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
  98       IF (EST.GT.H) GO TO 30
  99
 100       CALL GUFLD(XYZT,F)
 101       AT     = A + SECXS(1)
 102       BT     = B + SECYS(1)
 103       CT     = C + SECZS(1)
 104 *
 105       SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
 106       SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
 107       SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
 108       AT     = A + SECXS(2)
 109       BT     = B + SECYS(2)
 110       CT     = C + SECZS(2)
 111       SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
 112       SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
 113       SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
 114       DXT    = H * (A + SECXS(3))
 115       DYT    = H * (B + SECYS(3))
 116       DZT    = H * (C + SECZS(3))
 117       XT     = X + DXT
 118       YT     = Y + DYT
 119       ZT     = Z + DZT
 120       AT     = A + TWO*SECXS(3)
 121       BT     = B + TWO*SECYS(3)
 122       CT     = C + TWO*SECZS(3)
 123 *
 124       EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
 125       IF (EST.GT.2.*ABS(H)) GO TO 30
 126
 127       CALL GUFLD(XYZT,F)
 128 *
 129       Z      = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
 130       Y      = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
 131       X      = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
 132 *
 133       SECXS(4) = (BT*F(3) - CT*F(2))* PH2
 134       SECYS(4) = (CT*F(1) - AT*F(3))* PH2
 135       SECZS(4) = (AT*F(2) - BT*F(1))* PH2
 136       A      = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
 137       B      = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
 138       C      = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
 139 *
 140       EST    = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
 141      ++        ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
 142      ++        ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
 143 *
 144       IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
 145       ITER = ITER + 1
 146       NCUT = 0
 147 *               If too many iterations, go to HELIX
 148       IF (ITER.GT.MAXIT) GO TO 40
 149 *
 150       TL = TL + H
 151       IF (EST.LT.(DLT32)) THEN
 152          H = H*TWO
 153       ENDIF
 154       CBA    = ONE/ SQRT(A*A + B*B + C*C)
 155       VOUT(1) = X
 156       VOUT(2) = Y
 157       VOUT(3) = Z
 158       VOUT(4) = CBA*A
 159       VOUT(5) = CBA*B
 160       VOUT(6) = CBA*C
 161       REST = STEP - TL
 162       IF (STEP.LT.0.) REST = -REST
 163       IF (REST .GT. 1.E-5*ABS(STEP)) GO TO 20
 164 *
 165       GO TO 999
 166 *
 167 **              CUT STEP
 168    30 NCUT = NCUT + 1
 169 *               If too many cuts , go to HELIX
 170       IF (NCUT.GT.MAXCUT)       GO TO 40
 171       H = H*HALF
 172       GO TO 20
 173 *
 174 **              ANGLE TOO BIG, USE HELIX
 175    40 F1  = F(1)
 176       F2  = F(2)
 177       F3  = F(3)
 178       F4  = SQRT(F1**2+F2**2+F3**2)
 179       RHO = -F4*PINV
 180       TET = RHO * STEP
 181       IF(TET.NE.0.) THEN
 182          HNORM = ONE/F4
 183          F1 = F1*HNORM
 184          F2 = F2*HNORM
 185          F3 = F3*HNORM
 186 *
 187          HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
 188          HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
 189          HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)
 190
 191          HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
 192 *
 193          RHO1 = ONE/RHO
 194          SINT = SIN(TET)
 195          COST = TWO*SIN(HALF*TET)**2
 196 *
 197          G1 = SINT*RHO1
 198          G2 = COST*RHO1
 199          G3 = (TET-SINT) * HP*RHO1
 200          G4 = -COST
 201          G5 = SINT
 202          G6 = COST * HP
 203
 204          VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
 205          VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
 206          VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
 207
 208          VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
 209          VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
 210          VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
 211 *
 212       ELSE
 213          VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
 214          VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
 215          VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
 216 *
 217       ENDIF
 218 *
 219   999 END