GEANT321/neutron/seclf1.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 SECLF1(FSE,IFSE,EX,U,E)
  13 C        THIS ROUTINE SAMPLES AN EXIT ENERGY FROM
  14 C        A TABULATED DISTRIBUTION
  15 #include "geant321/minput.inc"
  16       DIMENSION FSE(*),IFSE(*)
  17       SAVE
  18       C=0.
  19       IP=1
  20       NRE=IFSE(IP)
  21       NE=IFSE(IP+1)
  22       IP=IP+2*NRE+1
  23 C       FIND THE TWO INCIDENT ENERGY DISTRIBUTIONS THAT BOUND E
  24 C       INCIDENT ENERGY IS BELOW THE FIRST INCIDENT ENERGY GIVEN
  25 C       USE THE FIRST DISTRIBUTION
  26       IP=IP+1
  27       IE=1
  28       E1=FSE(IP)
  29       IP1=IP
  30       IF(E.GT.E1)GO TO 10
  31       IPR=IP+1
  32       NR=IFSE(IPR)
  33       NP=IFSE(IPR+1)
  34       GO TO 50
  35    10 IP=IP+1
  36       IPR1=IP
  37       NP1=IFSE(IP+1)
  38       IP=IP+2*IFSE(IPR1)+1
  39       IP=IP+2*NP1
  40    20 IE=IE+1
  41       IP=IP+1
  42 C       INCIDENT ENERGY IS ABOVE THE LAST INCIDENT ENERGY GIVEN
  43 C       USE THE LAST DISTRIBUTION
  44       IF(IE.GT.NE)GO TO 40
  45       E2=FSE(IP)
  46       IF(E.LE.E2)GO TO 30
  47       E1=E2
  48       IP1=IP
  49       IP=IP+1
  50       IPR1=IP
  51       NP1=IFSE(IP+1)
  52       IP=IP+2*IFSE(IPR1)+1
  53       IP=IP+2*NP1
  54       GO TO 20
  55    30 IP2=IP
  56       IP=IP+1
  57       IPR2=IP
  58       NP2=IFSE(IP+1)
  59       IP=IP+2*IFSE(IPR2)+1
  60 C       DETERMINE THE INTERPOLATING SCHEME
  61       CALL INTSCH(IFSE,IE,IS,NRE)
  62 C       SELECT THE DISTRIBUTION TO SAMPLE FROM
  63       R=FLTRNF(0)
  64 C       INTERPOLATION SCHEMES OF 1 (CONSTANT) OR 2 (LINEAR) ALLOWED
  65       IF(IS.GT.2)GO TO 110
  66       PROB=(E2-E)/(E2-E1)
  67       IF(IS.EQ.1)PROB=1.0
  68       IF(R.LE.PROB)GO TO 40
  69 C       SELECT FROM THE SECOND DISTRIBUTION
  70       NP=NP2
  71       IP=IP2
  72       IPR=IPR2
  73       GO TO 50
  74 C       SELECT FROM THE FIRST DISTRIBUTION
  75 C       OR FROM THE LAST INCIDENT ENERGY
  76    40 NP=NP1
  77       IP=IP1
  78       IPR=IPR1
  79 C       SELECT THE EXIT ENERGY FROM THE TABULATED DISTRIBUTION
  80    50 CONTINUE
  81       ITRY = 0
  82    60 CONTINUE
  83       PROB=0.
  84       R=FLTRNF(0)
  85       NR=2*IFSE(IPR)+1
  86       DO 90  I=1,NP
  87          CALL INTSCH(IFSE(IPR),NP,IS,IFSE(IPR))
  88          N=IP+NR+1+2*I
  89          PROB1=PROB
  90          IF(I.EQ.1)GO TO 90
  91          IF(IS.EQ.1)GO TO 70
  92          IF(IS.GT.2)GO TO 110
  93          PROB=PROB+(FSE(N)+FSE(N-2))*(FSE(N-1)-FSE(N-3))/2.
  94          GO TO 80
  95    70    PROB=PROB+FSE(N-2)*(FSE(N-1)-FSE(N-3))
  96    80    IF(R.LE.PROB)GO TO 100
  97    90 CONTINUE
  98       ITRY = ITRY + 1
  99       IF(ITRY.LT.5) GOTO 60
 100       IF(R.LT..998)GO TO 120
 101   100 EX=FSE(N-3)+(R-PROB1)*(FSE(N-1)-FSE(N-3))/(PROB-PROB1)
 102       RETURN
 103   110 WRITE(IOUT,10000)IS
 104 10000 FORMAT(' MICAP: INTERPOLATION SCHEME=',I3,' IN SECLF1')
 105       GOTO 130
 106   120 WRITE(IOUT,10100)R,PROB
 107 10100 FORMAT(' MICAP: EXIT ENERGY NOT SELECTED IN SECLF1',1P2E13.5)
 108   130 WRITE(6,*) ' CALOR: ERROR in SECLF1 =====> STOP '
 109       STOP
 110       END