SUBROUTINE MUDECY
 
C-----------------------------------------------------------------------
C  MU(ON) DEC(A)Y
C
C  TREATES DECAY OF MUON INTO ELECTRON (INCLUDING POLARISATION EFFECTS)
C  INCLUDING NEUTRINOS, IF SELECTED
C  THIS SUBROUTINE IS CALLED FROM MUTRAC
C
C  DESIGN  : D. HECK    IK3  FZK KARLSRUHE
C-----------------------------------------------------------------------
 
      IMPLICIT NONE
*KEEP,CONST.
      COMMON /CONST/   PI,PI2,OB3,TB3,ENEPER
      DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
*KEEP,GENER.
      COMMON /GENER/   GEN,ALEVEL
      DOUBLE PRECISION GEN,ALEVEL
*KEEP,PAM.
      COMMON /PAM/     PAMA,SIGNUM
      DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
*KEEP,PARPAR.
      COMMON /PARPAR/  CURPAR,SECPAR,PRMPAR,OUTPAR,C,
     *                 E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
      DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
     *                 C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
      INTEGER          ITYPE,LEVL
*KEEP,PARPAE.
      DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
      EQUIVALENCE      (CURPAR(2),GAMMA),  (CURPAR(3),COSTHE),
     *                 (CURPAR(4), PHI ),  (CURPAR(5), H    ),
     *                 (CURPAR(6), T   ),  (CURPAR(7), X    ),
     *                 (CURPAR(8), Y   ),  (CURPAR(9), CHI  ),
     *                 (CURPAR(10),BETA),  (CURPAR(11),GCM  ),
     *                 (CURPAR(12),ECM )
*KEEP,POLAR.
      COMMON /POLAR/   POLART,POLARF
      DOUBLE PRECISION POLART,POLARF
*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
*KEND.
 
      DOUBLE PRECISION AUX2,COSDE,COSTH3,COS3CM,COS3C1,COS3C2,
     *                 E3CM,GAMMA3,PHI3CM,PHI3C2,PHI31,
     *                 P3CM,XI
      INTEGER          I
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
  444 FORMAT(' MUDECY: CURPAR=',1P,9E10.3)
 
C  MUON DECAYS INTO ELECTRON AND NEUTRINOS
      XI     = 2*ITYPE - 11
C  ELECTRON ENERGY SPECTRUM  N(E) * DE = CONST * E**2 * (3/2*E0-E) * DE
C  IS GAINED BY THE REJECTION/REFLECTION METHOD
   6  CALL RMMAR( RD,4,1 )
      IF ( RD(1)**2*(3.-RD(1)*2.)  .LT.  RD(2) )  RD(1) = 1.-RD(1)
      E3CM   = PAMA(2) + RD(1) * ( C(8) - PAMA(2) )
      IF ( E3CM .GT. 0.5D0*PAMA(5) )  GOTO 6
      P3CM   = SQRT( E3CM**2 - PAMA(2)**2 )
C  NOW DETERMINE COS3C1 AND PHI31 BY RANDOM SELECTION
C  WITH RESPECT TO THE POLARIZATION DIRECTION OF THE MUON IN THE MU CM
C  GIVEN BY POLART, POLARF
      COSDE  = 2.D0 * RD(4) - 1.D0
      AUX2   = ( 1. - 2.*RD(1) ) / ( 3. - 2.*RD(1) )
      IF ( ABS(AUX2) .GT. 1.D-2 ) THEN
        COS3C1 = XI*(SQRT(1.D0-(2.D0*COSDE-AUX2)*AUX2) - 1.D0) / AUX2
      ELSE
        COS3C1 = -XI * COSDE
      ENDIF
      PHI31  = RD(3)*PI2
 
C  NOW ADD ELECTRON EMISSION ANGLE COS3C1 TO THE POLARISATION DIRECTION
C  TO GET THE DIRECTION (RELATIVE TO THE CORSIKA COORDINATE SYSTEM)
      CALL ADDANG( POLART,POLARF, COS3C1,PHI31, COS3C2,PHI3C2 )
C  GET THE ELECTRON DIRECTION RELATIVE TO THE MUON LAB DIRECTION
      CALL ADDANI( CURPAR(3),CURPAR(4), COS3C2,PHI3C2, COS3CM,PHI3CM )
C  LORENTZ TRANSFORMATION TO THE LAB SYSTEM
      GAMMA3 = GAMMA * ( E3CM + BETA * P3CM * COS3CM ) / PAMA(2)
      COSTH3 = MIN( 1.D0, GAMMA * (P3CM * COS3CM + BETA * E3CM) /
     *                             (PAMA(2) * SQRT(GAMMA3**2 - 1.D0)) )
      CALL ADDANG( CURPAR(3),CURPAR(4), COSTH3,PHI3CM,
     *                                    SECPAR(3),SECPAR(4) )
      IF ( SECPAR(3) .GE. C(29) ) THEN
        SECPAR(1) = ITYPE - 3
        SECPAR(2) = GAMMA3
        DO 10  I = 5,8
          SECPAR(I) = CURPAR(I)
  10    CONTINUE
        SECPAR( 9)  = GEN
        SECPAR(10)  = ALEVEL
        CALL TSTACK
      ENDIF
      POLART = 0.D0
      POLARF = 0.D0
 
      RETURN
      END