source: trunk/MagicSoft/Simulation/Corsika/Mmcs/compt.f@ 18679

      SUBROUTINE COMPT
C                                VERSION 4.00  --  26 JAN 1986/1900
C******************************************************************
C   BUTCHER AND MESSEL'S CROSS SECTION EXPRESSION IS USED
C   (BUTCHER AND MESSEL, OP.CIT., P. 17-19,25), BUT THE
C   1/EPSILON PART IS NOT SAMPLED IN THE WAY THAT THEY DO.
C   THIS ROUTINE CALLS THEIR 'EPSILON' VARIABLE BY THE NAME 'BR'.
C   BR=FINAL PHOTON ENERGY /INITIAL PHOTON ENERGY.
C   BR0 = MINIMUM BR = 1./(1.+2.*(E(NP)/RM))
C   MAXIMUM BR IS 1.
C   BUTCHER AND MESSEL'S EXPRESSION FOR THE DIFFERENTIAL CROSS
C   SECTION IS PROPORTIONAL TO
C        (1./BR+BR)*(1.-BR*SINTHE**2/(1.+BR*BR))
C   WE SHALL SAMPLE FROM THE FIRST FACTOR FROM THE INTERVAL (BR0,1)
C   AND USE THE SECOND FACTOR AS A REJECTION FUNCTION.
C******************************************************************
      DOUBLE PRECISION PEIG,PESG,PESE
*KEEP,RANDPA.
      COMMON /RANDPA/  FAC,U1,U2,RD,NSEQ,ISEED,KNOR
      DOUBLE PRECISION FAC,U1,U2
      REAL             RD(3000)
      INTEGER          ISEED(103,10),NSEQ
      LOGICAL          KNOR
*KEEP,RUNPAR.
      COMMON /RUNPAR/  FIXHEI,THICK0,HILOECM,HILOELB,
     *                 STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
     *                 MONIOU,MDEBUG,NUCNUC,
     *                 CETAPE,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE,
     *                 DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
      COMMON /RUNPAC/  DSN,HOST,USER
      DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
      REAL             STEPFC
      INTEGER          NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE
      INTEGER          CETAPE
      CHARACTER*79     DSN
      CHARACTER*20     HOST,USER
 
      LOGICAL          DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
*KEEP,STACKE.
      COMMON/STACKE/   E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
      DOUBLE PRECISION E(60),TIME(60)
      REAL             X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
      INTEGER          IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
*KEND.
      COMMON/THRESH/RMT2,RMSQ,ESCD2,AP,API,AE,UP,UE,TE,THMOLL
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC
      COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL
      COMMON/ACLOCK/NCLOCK,JCLOCK
C_____IF (NCLOCK.GT.JCLOCK) THEN
C______WRITE(MDEBUG,* )' COMPT: NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
C______CALL AUSGB2
C_____END IF
      PEIG=E(NP)
      EIG=PEIG
      EGP=EIG*RMI
      BR0I=1.+2.*EGP
      BR0=1./BR0I
      ALPH1=LOG(BR0I)
      ALPH2=EGP*(BR0I+1.)*BR0*BR0
      SUMALP = ALPH1+ALPH2
371   CONTINUE
       CALL RMMAR(RNNO15,1,2)
       IF (ALPH1.GE.SUMALP*RNNO15) THEN
        CALL RMMAR(RNNO16,1,2)
        BR=EXP(ALPH1*RNNO16)*BR0
       ELSE
        CALL RMMAR(RD,2,2)
        BRP=RD(1)
        RNNO18=RD(2)
        IF (EGP.GE.(EGP+1.)*RNNO18) THEN
         CALL RMMAR(RNNO19,1,2)
         BRP=MAX(BRP,RNNO19)
        END IF
        BR=((BR0I-1.)*BRP+1.)*BR0
       END IF
       ESG=BR*EIG
       A1MIBR = 1.-BR
       TEMP=RM*A1MIBR/ESG
       SINTHE=MAX(0.0,TEMP*(2.0-TEMP))
       CALL RMMAR(RNNO20,1,2)
       IF(((1.0-RNNO20)*(1.0+BR*BR).GE.BR*SINTHE))GO TO372
      GO TO 371
372   CONTINUE
      SINTHE=SQRT(SINTHE)
      COSTHE=1.0-TEMP
      PESG=ESG
      PESE=PEIG-PESG+PRM
      ESE=PESE
      CALL UPHI(2,1)
      NP=NP+1
      PSQ=ESE*ESE-RMSQ
      IF (PSQ.LE.0.0) THEN
       COSTHE=0.
       SINTHE=-1.
      ELSE
       COSTHE=(ESE+ESG)*A1MIBR/SQRT(PSQ)
       SINTHE=-SQRT(MAX(0.0,1.0-COSTHE*COSTHE))
      END IF
      CALL UPHI(3,2)
      IF (ESE.LE.ESG) THEN
       IQ(NP)=3
       E(NP)=PESE
       E(NP-1)=PESG
      ELSE
       IQ(NP)=1
       IQ(NP-1)=3
       E(NP)=PESG
       E(NP-1)=PESE
       T=U(NP)
       U(NP)=U(NP-1)
       U(NP-1)=T
       T=V(NP)
       V(NP)=V(NP-1)
       V(NP-1)=T
       T=W(NP)
       W(NP)=W(NP-1)
       W(NP-1)=T
      END IF
      RETURN
      END