source: trunk/MagicSoft/Simulation/Corsika/Mmcs/mucoul.f

      SUBROUTINE MUCOUL(OMEGA,DENS,VSCAT)
 
C-----------------------------------------------------------------------
C   MU(ON) COUL(OMB SCATTERING OF SINGLE SCATTERING EVENTS)
C
C  TREATES SINGLE COULOMB SCATTERING FOR MUONS IN SMALL ANGLE
C  APPROXIMATION.
C  THIS SUBROUTINE IS IN ANALOGY WITH SUBROUTINE GMCOUL
C  (AUTHOR: G. LYNCH, LBL) OF GEANT321
C  SEE CERN PROGRAM LIBRARY LONG WRITEUP W5013
C  THIS SUBROUTINE IS CALLED FROM UPDATE
C  ARGUMENTS:
C   OMEGA = NUMBER OF SCATTERINGS FOR THE STEP
C   DENS  = LOCAL DENSITY
C   VSCAT = SCATTERING ANGLE
C
C  REDESIGN: 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,MUMULT.
      COMMON /MUMULT/  CHC,OMC,FMOLI
      DOUBLE PRECISION CHC,OMC
      LOGICAL          FMOLI
*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,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 DENS,OMCF,OMEGA,OMEGA0,PHIS,SUMX,SUMY,
     *                 THET,THMIN2,VSCAT
      INTEGER          I,NSCMX,NSCA
      DATA             OMCF/1.167D0/,NSCMX/50/
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,*)'MUCOUL: OMEGA=',SNGL(OMEGA),
     *                                    ' DENS=',SNGL(DENS)
 
C  COMPUTE NUMBER OF SCATTERS (POISSON DISTR. WITH MEAN OMEGA0)
      OMEGA0 = OMCF*OMEGA
      CALL MPOISS (OMEGA0,NSCA)
      IF ( NSCA .LE. 0 ) THEN
        VSCAT = 0.D0
        RETURN
      ENDIF
      NSCA = MIN(NSCA,NSCMX)
      CALL RMMAR(RD,2*NSCA,1)
 
C  THMIN2 IS THE SCREENING ANGLE
      THMIN2 = CHC**2/( OMCF*OMC * (PAMA(5)*BETA*GAMMA)**2 )
 
      SUMX = 0.D0
      SUMY = 0.D0
      DO 12 I = 1,NSCA
        THET  = SQRT( THMIN2*((1./RD(I)) - 1.) )
        PHIS  = PI2 * RD(NSCA+I)
        SUMX  = SUMX + THET*COS(PHIS)
        SUMY  = SUMY + THET*SIN(PHIS)
 12   CONTINUE
      VSCAT  = SQRT(SUMX**2 + SUMY**2)
 
      RETURN
      END