source: trunk/MagicSoft/Simulation/Corsika/Mmcs/prange.f@ 19094

      SUBROUTINE PRANGE(ARG)
 
C-----------------------------------------------------------------------
C  (DECAYING) P(ARTICLE'S) RANGE
C
C  DETERMINES MEAN FREE PATH FOR DECAYING PARTICLES
C  INCLUDING IONIZATION ENERGY LOSS,
C  FOR EACH LAYER OF THE ATMOSOHERE SEPARATELY
C  PRECISELY
C  THIS SUBROUTINE IS CALLED FROM BOX2
C  ARGUMENT:
C   ARG    = -LOG(RANDOM NUMBER) * SPEED OF LIGHT * LIFETIME
C
C  DESIGN  : D. HECK    IK3  FZK KARLSRUHE
C-----------------------------------------------------------------------
 
      IMPLICIT NONE
*KEEP,AIR.
      COMMON /AIR/     COMPOS,PROBTA,AVERAW,AVOGAD
      DOUBLE PRECISION COMPOS(3),PROBTA(3),AVERAW,AVOGAD
*KEEP,ATMOS.
      COMMON /ATMOS/   AATM,BATM,CATM,DATM
      DOUBLE PRECISION AATM(5),BATM(5),CATM(5),DATM(5)
*KEEP,ATMOS2.
      COMMON /ATMOS2/  HLAY,THICKL
      DOUBLE PRECISION HLAY(5),THICKL(5)
*KEEP,CONST.
      COMMON /CONST/   PI,PI2,OB3,TB3,ENEPER
      DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
*KEEP,OBSPAR.
      COMMON /OBSPAR/  OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
     *                 THETPR,PHIPR,NOBSLV
      DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
     *                 THETAP,THETPR(2),PHIP,PHIPR(2)
      INTEGER          NOBSLV
*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,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 AK,ARG,ARG0,BK,CHIT,DK,ELOSS
      DOUBLE PRECISION GAMK,GAMNEW,GAMSQ,GAM0,GMSQM1,H0,TH0
      INTEGER          ILAY
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,444) ARG,THICKH
  444 FORMAT(' PRANGE: -LOG(RD)*C*TAU = ',1P,E10.3,' THICKH=',E10.3)
 
C  LOOK WITHIN WHICH LAYER THE PARTICLE STARTS
      IF     ( H .LE. HLAY(2) ) THEN
        ILAY = 1
        TH0  = THICKH
      ELSEIF ( H .LE. HLAY(3) ) THEN
        ILAY = 2
        TH0  = THICKH
      ELSEIF ( H .LE. HLAY(4) ) THEN
        ILAY = 3
        TH0  = THICKH
      ELSE
        ILAY = 4
        TH0  = MAX( THICKH, 2.D-4 )
      ENDIF
C  SET START VALUES FOR ITERATION
      ARG0 = ARG
      CHIT = 0.D0
      GAM0 = GAMMA
      H0   = H
 
  2   CONTINUE
      GAM0   = MAX( GAM0, 1.0001D0 )
      GAMSQ  = GAM0**2
      GMSQM1 = GAMSQ - 1.D0
C  ENERGY LOSS BY IONIZATION
      ELOSS  = SIGNUM(ITYPE)**2 * C(22) *
     *           ( GAMSQ * (LOG(GMSQM1) + C(23)) / GMSQM1 - 1.D0 )
      ELOSS  = ELOSS / (PAMA(ITYPE) * COSTHE )
      BK     = ELOSS * (TH0 - AATM(ILAY))
      DK     = GAM0 + BK
      AK     = ARG0 * DK * COSTHE * DATM(ILAY)
      IF ( AK .LT. 174.D0 ) THEN
C  LIMIT FOR EXPONENT (ON IBM COMPUTER)
        GAMNEW = MAX( GAM0 * DK / ( GAM0 + BK * EXP(AK) ), 1.D0 )
      ELSE
        GAMNEW = 1.D0
      ENDIF
      GAMK   = GAM0 - ELOSS * ( THICKL(ILAY) - TH0)
      IF ( DEBUG ) WRITE(MDEBUG,*) 'PRANGE: GAMNEW,GAMK=',
     *                           SNGL(GAMNEW),SNGL(GAMK)
C  LOOK WETHER PARTICLE PENETRATES LAYER BOUNDARY OR DECAYS BEFORE
      IF ( GAMNEW .LT. GAMK  .AND.  ILAY. GT. 1 ) THEN
C  CALCULATE PORTION OF RANGE AND NEW START VALUES AT LAYER BOUNDARY
        ARG0 = ARG0 - ( H0 - HLAY(ILAY) + CATM(ILAY) * LOG(GAM0/GAMK) )
     *                       / (DK * COSTHE)
        CHIT = CHIT + (THICKL(ILAY) - TH0) / COSTHE
        GAM0 = GAMK
        H0   = HLAY(ILAY)
        TH0  = THICKL(ILAY)
        ILAY = ILAY - 1
        GOTO 2
      ENDIF
C  PENETRATED MATTER THICKNESS
      CHI = CHIT + (GAM0 - GAMNEW) / (ELOSS*COSTHE)
      IF ( DEBUG ) WRITE(MDEBUG,445) CHI
  445 FORMAT(' PRANGE: CHI = ',1P,E10.3)
 
      RETURN
      END