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

      SUBROUTINE MSCAT
C                                VERSION 4.00  --  26 JAN 1986/1900
C******************************************************************
      COMMON/ELECIN/EKELIM,ICOMP,EKE0,EKE1,CMFP0,CMFP1,RANGE0,RANGE1, XR
     *0,TEFF0,BLCC,XCC,PICMP0(1),PICMP1(1),EICMP0(1),EICMP1(1),MPEEM(1),
     * ESIG0(500),ESIG1(500),PSIG0(500),PSIG1(500),EDEDX0(500),EDEDX1(50
     *0),PDEDX0(500),PDEDX1(500),EBR10(500),EBR11(500),PBR10(500),PBR11(
     *500),PBR20(500),PBR21(500),TMXS0(500),TMXS1(500),CMFPE0(1),CMFPE1(
     *1),CMFPP0(1),CMFPP1(1),ERANG0(1),ERANG1(1),PRANG0(1),PRANG1(1),CXC
     *2E0(1),CXC2E1(1),CXC2P0(1),CXC2P1(1),CLXAE0(1),CLXAE1(1),CLXAP0(1)
     *,CLXAP1(1), THR0(1,1),THR1(1,1),THR2(1,1),THRI0(1,1),THRI1(1,1),TH
     *RI2(1,1),FSTEP(16),FSQR(16),MSMAP(200), VERT1(1000),VERT2(100,16),
     *MSTEPS,JRMAX,MXV1, MXV2,NBLC,NRNTH,NRNTHI,BLC0,BLC1,RTHR0,RTHR1,RT
     *HRI0,RTHRI1
*KEEP,EPCONT.
      COMMON/EPCONT/   EDEP,RATIO,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,IDISC,
     *                 IROLD,IRNEW,RHOFAC, EOLD,ENEW,EKE,ELKE,BETA2,GLE,
     *                 TSCAT,IAUSFL
      DOUBLE PRECISION EDEP,RATIO
      REAL             TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,RHOFAC,EOLD,ENEW,
     *                 EKE,ELKE,BETA2,GLE,TSCAT
      INTEGER          IDISC,IROLD,IRNEW,IAUSFL(29)
*KEND.
      COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(6),RHOR(6),IRAYLR(6)
      COMMON/MULTS/NG21,B0G21,B1G21,G210(7),G211(7),G212(7), NG22,B0G22,
     *B1G22,G220(8),G221(8),G222(8), NG31,B0G31,B1G31,G310(11),G311(11),
     *G312(11), NG32,B0G32,B1G32,G320(25),G321(25),G322(25), NBGB,B0BGB,
     *B1BGB,BGB0(8),BGB1(8),BGB2(8)
*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/UPHIIN/SINC0,SINC1,SIN0(20002),SIN1(20002)
      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,* )' MSCAT: NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
C______CALL AUSGB2
C_____END IF
      VSTEFF=TVSTEP*RHOFAC
      OMEGA0=BLCC*VSTEFF/BETA2
      IF ((OMEGA0.LE.1.0)) THEN
       SINTHE=0.0
       COSTHE=1.0
       THETA=0.0
       NOSCAT=NOSCAT+1
       RETURN
      END IF
      BLC=LOG(OMEGA0)
      IF (BLC.LE.1.30685) THEN
       B=1.530394*BLC
      ELSE
       IB=B0BGB+BLC*B1BGB
       IF (IB.GT.NBGB) THEN
        WRITE(KMPO,940)IB
940     FORMAT('MSCAT: NBGB<IB=',I5)
       END IF
       B=BGB0(IB)+BLC*(BGB1(IB)+BLC*BGB2(IB))
      END IF
      XR=XCC*SQRT(MAX(0.,VSTEFF*B))/(EOLD*BETA2)
      IF (B.GT.2.0) THEN
       BI=1./B
       BMD=1./(1.0+1.75*BI)
       BM1=(1.0-2.0*BI)*BMD
       BM2=(1.0+0.025*BI)*BMD
      ELSE
       BI=0.5
       BM1=(1.-2./B)*0.53333333
       BM2=0.54
      END IF
951   CONTINUE
       CALL RMMAR(RMS1,1,2)
       IF (RMS1.LE.BM1) THEN
        CALL RMMAR(RMS2,1,2)
        IF((RMS2.EQ.0.0))RMS2=1.E-30
        THR=SQRT(MAX(0.,-LOG(RMS2)))
       ELSE IF((RMS1.LE.BM2)) THEN
        CALL RMMAR(RD,3,2)
        RMS3=RD(1)
        RMS4=RD(2)
        RMS5=RD(3)
        ETA=MAX(RMS3,RMS4)
        I31=B0G31+ETA*B1G31
        G31=G310(I31)+ETA*(G311(I31)+ETA*G312(I31))
        I32=B0G32+ETA*B1G32
        G32=G320(I32)+ETA*(G321(I32)+ETA*G322(I32))
        G3=G31+G32*BI
        IF((RMS5.GT.G3))GO TO951
        THR=1.0/ETA
       ELSE
        CALL RMMAR(RD,2,2)
        RMS6=RD(1)
        RMS7=RD(2)
        THR=RMS6
        I21=B0G21+THR*B1G21
        G21=G210(I21)+THR*(G211(I21)+THR*G212(I21))
        I22=B0G22+THR*B1G22
        G22=G220(I22)+THR*(G221(I22)+THR*G222(I22))
        G2=G21+G22*BI
        IF((RMS7.GT.G2))GO TO951
       END IF
       THETA=THR*XR
       IF((THETA.GE.PI))GO TO951
       LTHETA=SINC1*THETA+SINC0
       SINTHE=SIN1(LTHETA)*THETA+SIN0(LTHETA)
       CALL RMMAR(RMS8,1,2)
       IF(((RMS8*RMS8*THETA.LE.SINTHE)))GO TO952
      GO TO 951
952   CONTINUE
      CTHET=PI5D2-THETA
      LCTHET=SINC1*CTHET+SINC0
      COSTHE=SIN1(LCTHET)*CTHET+SIN0(LCTHET)
      RETURN
      END