GEANT321/gphys/gcoeff.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:23  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GCOEFF
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *  Calculates the coefficients for the energy loss               *
  17 C.    *     interpolation                                              *
  18 C.    *     There are 4 tables : electron,positron,muon,proton         *
  19 C.    *                                                                *
  20 C.    *    ==>Called by : GPHYSI                                       *
  21 C.    *       Author      F.Carminati *********                        *
  22 C.    *                                                                *
  23 C.    ******************************************************************
  24 C.
  25 #include "geant321/gcbank.inc"
  26 #include "geant321/gctrak.inc"
  27 #include "geant321/gcjloc.inc"
  28 #include "geant321/gcmulo.inc"
  29 #include "geant321/gconsp.inc"
  30 #include "geant321/gcmate.inc"
  31 #if !defined(CERNLIB_SINGLE)
  32       DOUBLE PRECISION CX1,CX2,CX3,CY1,CY2,CY3,CDEN1,CDEN2,CDEN3
  33       DOUBLE PRECISION ACOEFF,BCOEFF,CCOEFF,XRAT,CCOEF1,CCOEF3
  34       DOUBLE PRECISION SQEPSM,CFACT
  35       PARAMETER (EPSMAC=1E-6)
  36 #endif
  37 #if defined(CERNLIB_SINGLE)
  38       PARAMETER (EPSMAC=1E-11)
  39 #endif
  40 *
  41       SQEPSM = MAX(1.,91./NEK1)*10.*SQRT(EPSMAC)
  42       DO 10 IEKBIN=1,NEK1-2
  43 *
  44          I1 = IEKBIN
  45          I2 = I1 + 1
  46          I3 = I2 + 1
  47          CY1 = ELOW(I1)
  48          CY2 = ELOW(I2)
  49          CY3 = ELOW(I3)
  50          IECOEF = 3*(IEKBIN-1)
  51 *
  52 * *** Electrons
  53 *
  54          JRANG = LQ(JMA-15)
  55          JCOEF = LQ(JMA-17)
  56 *
  57          CX1 = Q(JRANG+I1)
  58          CX2 = Q(JRANG+I2)
  59          CX3 = Q(JRANG+I3)
  60          IF(CX1.NE.CX2.AND.CX1.NE.CX3.AND.CX2.NE.CX3) THEN
  61             CDEN1 = 1./((CX1-CX2)*(CX1-CX3))
  62             CDEN2 = 1./((CX2-CX1)*(CX2-CX3))
  63             CDEN3 = 1./((CX3-CX1)*(CX3-CX2))
  64             ACOEFF = CY1*CDEN1+CY2*CDEN2+CY3*CDEN3
  65             BCOEFF = -(CY1*(CX2+CX3)*CDEN1+CY2*(CX1+CX3)*CDEN2+
  66      +                 CY3*(CX1+CX2)*CDEN3)
  67             CCOEFF = CY1*CX2*CX3*CDEN1+CX1*CY2*CX3*CDEN2+
  68      +               CX1*CX2*CY3*CDEN3
  69             IF(ACOEFF.EQ.0.) THEN
  70                XRAT = 0.
  71             ELSEIF(BCOEFF.GT.0.) THEN
  72                CFACT  = SQRT(ABS(ACOEFF))
  73                CCOEF1 = SQRT(ABS(CCOEFF-CY1))*CFACT
  74                CCOEF3 = SQRT(ABS(CCOEFF-CY3))*CFACT
  75                XRAT   = MAX(CCOEF1,CCOEF3)/BCOEFF
  76             ELSE
  77                XRAT=1.
  78             ENDIF
  79             IF(XRAT.LE.SQEPSM) THEN
  80                Q(JCOEF+IECOEF+1) = 0.
  81                Q(JCOEF+IECOEF+2) = BCOEFF
  82                Q(JCOEF+IECOEF+3) = CCOEFF
  83             ELSE
  84                Q(JCOEF+IECOEF+1) = ACOEFF
  85                Q(JCOEF+IECOEF+2) = 0.5*BCOEFF/ACOEFF
  86                Q(JCOEF+IECOEF+3) = CCOEFF/ACOEFF
  87             ENDIF
  88          ENDIF
  89 *
  90 * *** Positons
  91 *
  92          JRANG = LQ(JMA-15) + NEK1
  93          JCOEF = LQ(JMA-17) +3*NEK1
  94 *
  95          CX1 = Q(JRANG+I1)
  96          CX2 = Q(JRANG+I2)
  97          CX3 = Q(JRANG+I3)
  98          IF(CX1.NE.CX2.AND.CX1.NE.CX3.AND.CX2.NE.CX3) THEN
  99             CDEN1 = 1./((CX1-CX2)*(CX1-CX3))
 100             CDEN2 = 1./((CX2-CX1)*(CX2-CX3))
 101             CDEN3 = 1./((CX3-CX1)*(CX3-CX2))
 102             ACOEFF = CY1*CDEN1+CY2*CDEN2+CY3*CDEN3
 103             BCOEFF = -(CY1*(CX2+CX3)*CDEN1+CY2*(CX1+CX3)*CDEN2+
 104      +                 CY3*(CX1+CX2)*CDEN3)
 105             CCOEFF = CY1*CX2*CX3*CDEN1+CX1*CY2*CX3*CDEN2+
 106      +               CX1*CX2*CY3*CDEN3
 107             IF(ACOEFF.EQ.0.) THEN
 108                XRAT = 0.
 109             ELSEIF(BCOEFF.GT.0.) THEN
 110                CFACT  = SQRT(ABS(ACOEFF))
 111                CCOEF1 = SQRT(ABS(CCOEFF-CY1))*CFACT
 112                CCOEF3 = SQRT(ABS(CCOEFF-CY3))*CFACT
 113                XRAT   = MAX(CCOEF1,CCOEF3)/BCOEFF
 114             ELSE
 115                XRAT=1.
 116             ENDIF
 117             IF(XRAT.LE.SQEPSM) THEN
 118                Q(JCOEF+IECOEF+1) = 0.
 119                Q(JCOEF+IECOEF+2) = BCOEFF
 120                Q(JCOEF+IECOEF+3) = CCOEFF
 121             ELSE
 122                Q(JCOEF+IECOEF+1) = ACOEFF
 123                Q(JCOEF+IECOEF+2) = 0.5*BCOEFF/ACOEFF
 124                Q(JCOEF+IECOEF+3) = CCOEFF/ACOEFF
 125             ENDIF
 126          ENDIF
 127 *
 128 * *** Muons
 129 *
 130          JRANG = LQ(JMA-16)
 131          JCOEF = LQ(JMA-18)
 132 *
 133          CX1 = Q(JRANG+I1)
 134          CX2 = Q(JRANG+I2)
 135          CX3 = Q(JRANG+I3)
 136          IF(CX1.NE.CX2.AND.CX1.NE.CX3.AND.CX2.NE.CX3) THEN
 137             CDEN1 = 1./((CX1-CX2)*(CX1-CX3))
 138             CDEN2 = 1./((CX2-CX1)*(CX2-CX3))
 139             CDEN3 = 1./((CX3-CX1)*(CX3-CX2))
 140             ACOEFF = CY1*CDEN1+CY2*CDEN2+CY3*CDEN3
 141             BCOEFF = -(CY1*(CX2+CX3)*CDEN1+CY2*(CX1+CX3)*CDEN2+
 142      +                 CY3*(CX1+CX2)*CDEN3)
 143             CCOEFF = CY1*CX2*CX3*CDEN1+CX1*CY2*CX3*CDEN2+
 144      +               CX1*CX2*CY3*CDEN3
 145             IF(ACOEFF.EQ.0.) THEN
 146                XRAT = 0.
 147             ELSEIF(BCOEFF.GT.0.) THEN
 148                CFACT  = SQRT(ABS(ACOEFF))
 149                CCOEF1 = SQRT(ABS(CCOEFF-CY1))*CFACT
 150                CCOEF3 = SQRT(ABS(CCOEFF-CY3))*CFACT
 151                XRAT   = MAX(CCOEF1,CCOEF3)/BCOEFF
 152             ELSE
 153                XRAT=1.
 154             ENDIF
 155             IF(XRAT.LE.SQEPSM) THEN
 156                Q(JCOEF+IECOEF+1) = 0.
 157                Q(JCOEF+IECOEF+2) = BCOEFF
 158                Q(JCOEF+IECOEF+3) = CCOEFF
 159             ELSE
 160                Q(JCOEF+IECOEF+1) = ACOEFF
 161                Q(JCOEF+IECOEF+2) = 0.5*BCOEFF/ACOEFF
 162                Q(JCOEF+IECOEF+3) = CCOEFF/ACOEFF
 163             ENDIF
 164          ENDIF
 165 *
 166 * *** Protons
 167 *
 168          JRANG = LQ(JMA-16) + NEK1
 169          JCOEF = LQ(JMA-18) +3*NEK1
 170 *
 171          CX1 = Q(JRANG+I1)
 172          CX2 = Q(JRANG+I2)
 173          CX3 = Q(JRANG+I3)
 174          IF(CX1.NE.CX2.AND.CX1.NE.CX3.AND.CX2.NE.CX3) THEN
 175             CDEN1 = 1./((CX1-CX2)*(CX1-CX3))
 176             CDEN2 = 1./((CX2-CX1)*(CX2-CX3))
 177             CDEN3 = 1./((CX3-CX1)*(CX3-CX2))
 178             ACOEFF = CY1*CDEN1+CY2*CDEN2+CY3*CDEN3
 179             BCOEFF = -(CY1*(CX2+CX3)*CDEN1+CY2*(CX1+CX3)*CDEN2+
 180      +                 CY3*(CX1+CX2)*CDEN3)
 181             CCOEFF = CY1*CX2*CX3*CDEN1+CX1*CY2*CX3*CDEN2+
 182      +               CX1*CX2*CY3*CDEN3
 183             IF(ACOEFF.EQ.0.) THEN
 184                XRAT = 0.
 185             ELSEIF(BCOEFF.GT.0.) THEN
 186                CFACT  = SQRT(ABS(ACOEFF))
 187                CCOEF1 = SQRT(ABS(CCOEFF-CY1))*CFACT
 188                CCOEF3 = SQRT(ABS(CCOEFF-CY3))*CFACT
 189                XRAT   = MAX(CCOEF1,CCOEF3)/BCOEFF
 190             ELSE
 191                XRAT=1.
 192             ENDIF
 193             IF(XRAT.LE.SQEPSM) THEN
 194                Q(JCOEF+IECOEF+1) = 0.
 195                Q(JCOEF+IECOEF+2) = BCOEFF
 196                Q(JCOEF+IECOEF+3) = CCOEFF
 197             ELSE
 198                Q(JCOEF+IECOEF+1) = ACOEFF
 199                Q(JCOEF+IECOEF+2) = 0.5*BCOEFF/ACOEFF
 200                Q(JCOEF+IECOEF+3) = CCOEFF/ACOEFF
 201             ENDIF
 202          ENDIF
 203 *
 204    10 CONTINUE
 205 *
 206       END