*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:21:59  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.48  by  S.Giani
*-- Author :
      SUBROUTINE THRSEL(NE,NP7,NB7,E,EOUT,FM,TDK,ETH,ALPHA,BETA,W,
     +                  PMUPS,PMDS,F,AWR,IIN,IPTMD,IPMDS,IOUT)
C      THIS ROUTINE SELECTS THE EXIT ENERGY AND SCATTERING ANGLE
C      FROM S(ALPHA,BETA) DATA TABLES
      DIMENSION ETH(*),ALPHA(*),BETA(*),ABETA2(2),PMDS(*),
     +   IPTMD(NE),W(NE),IPMDS(*),PMUPS(NB7,NE),F(NP7,NB7),AWR(*)
      SAVE
C
      AAWR=AWR(IIN)
      DO 10 IE=1,NE
         IF(E.LT.ETH(IE))GO TO 20
   10 CONTINUE
      INDX=NE
      GO TO 30
   20 INDX=IE
      IF(INDX.EQ.1)GO TO 30
      R=FLTRNF(0)
      DELE=(ETH(INDX)-E)/(ETH(INDX)-ETH(INDX-1))
      DELEN=DELE
      GO TO 40
   30 DELEN=1.0
      INDX=2
   40 PROB=DELEN*W(INDX-1)+(1.0-DELEN)*W(INDX)
      IF(R.LE.PROB)GO TO 120
C       NEUTRON DOWNSCATTERS
      R=FLTRNF(0)
      DO 90 I=1,2
         IP=IPTMD(INDX)
         NB=IPMDS(IP)
         DO 50 IB=1,NB
            IF(R.LE.PMDS(IP+NB+IB))GO TO 60
   50    CONTINUE
         WRITE(IOUT,10000)PMDS(IP+2*NB)
10000 FORMAT(' MICAP: CUMULATIVE DOWNSCATTER DIST. DOES NOT END ',
     +       'IN 1.0 IN THRSEL',E12.4)
         PRINT *,' CALOR: ERROR in MICAP ====> STOP '
         STOP
   60    IF(IB.EQ.1)GO TO 70
         DELE=(PMDS(IP+NB+IB)-R)/(PMDS(IP+NB+IB)-PMDS(IP+NB+IB-1))
         ABETA=DELE*(PMDS(IP+IB-1)-PMDS(IP+IB))+PMDS(IP+IB)
         GO TO 80
   70    ABETA=BETA(IB)
   80    ABETA2(I)=ABETA
         INDX=INDX-1
   90 CONTINUE
      ABETA=DELEN*ABETA2(2)+(1.0-DELEN)*ABETA2(1)
      EOUT=E-TDK*ABETA
      IF(EOUT.LT.1.0E-05)EOUT=1.0E-05
      DO 100 IB=1,NB7
         IF(ABETA.LE.BETA(IB))GO TO 110
  100 CONTINUE
      IB=NB7
  110 DELE=(ABETA-BETA(IB))/(BETA(IB-1)-BETA(IB))
      GO TO 180
C       NEUTRON UPSCATTERS
  120 R=FLTRNF(0)
      DO 170 I=1,2
         DO 130 IB=1,NB7
            IF(R.LE.PMUPS(IB,INDX))GO TO 140
  130    CONTINUE
         WRITE(IOUT,10100)PMUPS(NB7,INDX)
10100 FORMAT(' MICAP: CUMULATIVE UPSCATTER DIST. DOES NOT END ',
     +       'IN 1.0 IN THRSEL',E12.4)
         PRINT *,' CALOR: ERROR in MICAP ====> STOP '
         STOP
  140    IF(IB.EQ.1)GO TO 150
         DELE=(PMUPS(IB,INDX)-R)/(PMUPS(IB,INDX)-PMUPS(IB-1,INDX))
         ABETA=DELE*(BETA(IB-1)-BETA(IB))+BETA(IB)
         GO TO 160
  150    WRITE(IOUT,10200)PMUPS(1,INDX)
10200 FORMAT(' MICAP: CUMULATIVE UPSCATTER DIST. DOES NOT BEGIN ',
     +       'AT 0.0 IN THRSEL',E12.4)
         PRINT *,' CALOR: ERROR in MICAP ====> STOP '
         STOP
  160    ABETA2(I)=ABETA
         INDX=INDX-1
  170 CONTINUE
      ABETA=DELEN*ABETA2(2)+(1.0-DELEN)*ABETA2(1)
      EOUT=E+TDK*ABETA
C       SELECT ANGLE
  180 AMAX=(EOUT+E+2.0*SQRT(E*EOUT))/(AAWR*TDK)
      AMIN=(EOUT+E-2.0*SQRT(E*EOUT))/(AAWR*TDK)
      DO 190 IA=1,NP7
         IF(AMAX.LT.ALPHA(IA))GO TO 200
  190 CONTINUE
      IA=NP7
      DELA=0.0
      GO TO 210
  200 DELA=(ALPHA(IA)-AMAX)/(ALPHA(IA)-ALPHA(IA-1))
  210 F4=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
      F3=DELE*(F(IA-1,IB-1)-F(IA-1,IB))+F(IA-1,IB)
      F2=DELA*(F3-F4)+F4
      DO 220 IA=1,NP7
         IF(AMIN.LT.ALPHA(IA))GO TO 230
  220 CONTINUE
      IA=NP7
      DELA=0.0
      GO TO 240
  230 DELA=(ALPHA(IA)-AMIN)/(ALPHA(IA)-ALPHA(IA-1))
  240 F4=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
      F3=DELE*(F(IA-1,IB-1)-F(IA-1,IB))+F(IA-1,IB)
      F1=DELA*(F3-F4)+F4
      R=FLTRNF(0)
      F0=R*F2+(1.0-R)*F1
      F1=0.0
      DO 250 IA=1,NP7
         F2=DELE*(F(IA,IB-1)-F(IA,IB))+F(IA,IB)
         IF(F0.LE.F2)GO TO 260
         F1=F2
  250 CONTINUE
  260 IF(F1.EQ.F2)GO TO 270
      DELA=(F2-F0)/(F2-F1)
      GO TO 280
  270 ALP=ALPHA(IA)
      GO TO 290
  280 ALP=DELA*ALPHA(IA-1)+(1.0-DELA)*ALPHA(IA)
  290 FM=(E+EOUT-ALP*AAWR*TDK)/(2.0*SQRT(E*EOUT))
      IF(ABS(FM).LE.1.0)RETURN
      WRITE(IOUT,10300)FM,E,EOUT,R,IA,IB
10300 FORMAT(' MICAP: ERROR IN THRSEL, COSINE OF ANGLE >1'/,
     +' ',1P4E12.4,2I11)
      FM=2.0*FLTRNF(0)-1.0
      RETURN
      END