GEANT321/neutron/seclf5.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:59  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/04 23/02/95  14.46.01  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE SECLF5(FSE,IFSE,EX,U,E)
  13 C       THIS ROUTINE SAMPLES AN EXIT ENERGY FROM
  14 C       A GENERAL EVAPORATION SPECTRUM
  15 #include "geant321/minput.inc"
  16       DIMENSION FSE(*),IFSE(*)
  17       SAVE
  18 C       DETERMINE THETA
  19       IP=1
  20       NR=IFSE(IP)
  21       NE=IFSE(IP+1)
  22       IP=2*NR+IP
  23       EMAX=E-U
  24       R=FLTRNF(0)
  25       DO 10 I=1,NE
  26          IP=IP+2
  27          IF(E.LE.FSE(IP))GO TO 20
  28    10 CONTINUE
  29       GO TO 80
  30    20 IF(I.EQ.1)GO TO 90
  31 C       DETERMINE THE INTERPOLATING SCHEME
  32       CALL INTSCH(IFSE,I,IS,NR)
  33       E2=FSE(IP)
  34       E1=FSE(IP-2)
  35       CALL INTERP(E,THETA,E1,FSE(IP-1),E2,FSE(IP+1),IS)
  36       IP=IP+2+(NE-I)*2
  37 C       DETERMINE X
  38    30 NF=IFSE(IP+1)
  39       NR=IFSE(IP)
  40       IPR=IP
  41       IP=IP+1+2*NR
  42       PROB=0.
  43       DO 60 I=1,NF
  44          N=IP+2*I
  45          PROB1=PROB
  46          CALL INTSCH(IFSE(IPR),I,IS,NR)
  47          IF(I.EQ.1)GO TO 60
  48          IF(IS.EQ.1)GO TO 40
  49          IF(IS.GT.2)GO TO 100
  50          PROB=PROB+(FSE(N)+FSE(N-2))*(FSE(N-1)-FSE(N-3))/2.
  51          GO TO 50
  52    40    PROB=PROB+FSE(N-2)*(FSE(N-1)-FSE(N-3))
  53    50    CONTINUE
  54          IF(R.LE.PROB)GO TO 70
  55    60 CONTINUE
  56    70 X=FSE(N-3)+(R-PROB1)*(FSE(N-1)-FSE(N-3))/(PROB-PROB1)
  57 C       SELECT THE EXIT ENERGY FROM THE GENERAL EVAPORATION SPECTRUM
  58       EX=THETA*X
  59       IF(EX.LE.EMAX)RETURN
  60 #if defined(CERNLIB_MDEBUG)
  61       WRITE(IOUT,10000)EX,EMAX
  62 10000 FORMAT(' MICAP: WARNING-EX,EMAX=',1P2E13.5,' IN SECLF5')
  63 #endif
  64       EX=EMAX
  65       RETURN
  66 C       INCIDENT ENERGY IS ABOVE THE LAST INCIDENT ENERGY GIVEN
  67 C       USE THE LAST DISTRIBUTION
  68    80 THETA=FSE(IP+1)
  69       IP=IP+2
  70       GO TO 30
  71 C       INCIDENT ENERGY IS BELOW THE FIRST INCIDENT ENERGY GIVEN
  72 C       USE THE FIRST DISTRIBUTION
  73    90 THETA=FSE(IP+1)
  74       IP=IP+2*(NE-I)+2
  75       GO TO 30
  76   100 WRITE(IOUT,10100)IS
  77 10100 FORMAT(' MICAP: INTERPOLATION SCHEME=',I3,' IN SECLF5')
  78       WRITE(6,*) ' CALOR: ERROR in SECLF5 =====> STOP '
  79       STOP
  80       END