GEANT321/gstrag/gkokri.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:37  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.25  by  S.Giani
  11 *-- Author :
  12       FUNCTION GKOKRI(E,EMINEV,EMAXEV)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *     Input energy in eV, Sandia tables in keV                   *
  17 C.    *                                                                *
  18 C.    *                                                                *
  19 C.    ******************************************************************
  20 C.
  21 #include "geant321/gcbank.inc"
  22 #include "geant321/gcjloc.inc"
  23 #include "geant321/gconsp.inc"
  24 #include "geant321/gcmate.inc"
  25 #include "geant321/gc10ev.inc"
  26 #if !defined(CERNLIB_SINGLE)
  27       DOUBLE PRECISION EINV,ECUR,ONE,RES,ZERO,EBEG,EEND
  28       DOUBLE PRECISION ALPH,BETA,GAMM,E1,E2,GPSCIN,POLE
  29       DOUBLE PRECISION EMAX,EMIN,REST
  30 #endif
  31       PARAMETER (ONE=1,ZERO=0)
  32 C.
  33 C.    ------------------------------------------------------------------
  34 C.
  35       REST   = ZERO
  36       EMAX   = EMAXEV*1E-3
  37       EMIN   = EMINEV*1E-3
  38       ECUR   = E*1E-3
  39       EINV   = ONE/ECUR
  40 C Use Sandia data
  41       JPHXS  = LQ(JPHOT-1)
  42       NZ     = Q(JPHXS+1)
  43       JWEIGH = JPHXS+1+2*NZ
  44       DO 30 JZ=1,NZ
  45          RES    = ZERO
  46          EBEG   = EMIN
  47          WEIGHT = Q(JWEIGH+JZ)
  48          JPHFN  = LQ(JPHXS-JZ)
  49          IPOINT = JPHFN+1
  50          IMAX   = Q(IPOINT)
  51          IPOINT = IPOINT-4
  52          DO 10 I = 1,IMAX
  53             IPOINT = IPOINT+5
  54             EEND   = Q(IPOINT)
  55             IF(EEND.GT.EMIN) THEN
  56                E1     = MAX(EBEG,EMIN)
  57                E2     = MIN(EEND,EMAX)
  58                J      = IPOINT+1
  59                IF(ECUR.LE.E2.AND.ECUR.GE.E1) THEN
  60 *
  61 * *** The pole of the integration is in this interval
  62                   EPS1 = (ECUR-EBEG)/ECUR
  63                   EPS2 = (EEND-ECUR)/EEND
  64                   IF(EPS1.LT.EPS2) THEN
  65 *
  66 * *** First the pole and then a simple integration
  67                      ALPH = ONE-EPS1
  68                      E1 = ECUR/ALPH
  69                      E2 = EEND
  70                   ELSE
  71 *
  72 * *** First a simple integration and then the pole
  73                      ALPH = ONE-EPS2
  74                      E1 = EBEG
  75                      E2 = ECUR*ALPH
  76                   ENDIF
  77                   BETA = -LOG(ALPH)
  78                   GAMM = ONE/ALPH
  79                   POLE = EINV*(Q(J)*BETA+EINV*( Q(J+1)*(GAMM-ALPH)+
  80      +            EINV*( Q(J+2)*(0.5*GAMM**2+BETA-0.5*ALPH**2)+EINV*
  81      +            Q(J+3)*(GAMM**3/3+GAMM-ALPH-ALPH**3/3))))
  82                   RES = RES - POLE
  83                ENDIF
  84 *
  85 * *** This is a normal integration
  86                RES = RES + GPSCIN(E1,E2,ECUR,Q(J))
  87             ENDIF
  88             EBEG = EEND
  89             IF(EBEG.GE.EMAX) GO TO 20
  90    10    CONTINUE
  91    20    REST = REST+WEIGHT*RES
  92    30 CONTINUE
  93 C RES value is in cm**2/(g keV)
  94       GKOKRI = REST*1E-6*A/AVO
  95 C Now in Megabarns/eV
  96 C
  97       END