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

Last change on this file was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 5.6 KB
Line 
1 SUBROUTINE EGS4( EEIN )
2
3C---------------------------------------------------------------------
4C E(LECTRON) G(AMMA) S(HOWER)
5C
6C TREATES ELECTROMAGNETIC SUBSHOWER
7C THIS SUBROUTINE PACKAGE IS CALLED FROM EM
8C ARGUMENT:
9C EEIN = (R8) INCOMING PARTICLE ENERGY (GEV)
10C
11C DESIGN : D. HECK IK3 FZK KARLSRUHE
12C---------------------------------------------------------------------
13
14 DOUBLE PRECISION EEIN
15*KEEP,GENER.
16 COMMON /GENER/ GEN,ALEVEL
17 DOUBLE PRECISION GEN,ALEVEL
18*KEND.
19 COMMON/GEOM/ZALTIT,BOUND(6),NEWOBS,OBSLVL(10)
20*KEEP,LONGI.
21 COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
22 * NSTEP,LLONGI,FLGFIT
23 DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
24 * SPLONG(0:1040,9),THSTEP,THSTPI
25 INTEGER NSTEP
26 LOGICAL LLONGI,FLGFIT
27*KEND.
28 COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(6),RHOR(6),IRAYLR(6)
29*KEEP,OBSPAR.
30 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
31 * THETPR,PHIPR,NOBSLV
32 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
33 * THETAP,THETPR(2),PHIP,PHIPR(2)
34 INTEGER NOBSLV
35*KEEP,PARPAR.
36 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
37 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
38 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
39 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
40 INTEGER ITYPE,LEVL
41*KEND.
42 DOUBLE PRECISION PI0MSQ
43 COMMON/PION/PI0MSQ,PITHR,PICMAS,PI0MAS,AMASK0,AMASKC,AMASPR,AMASNT
44*KEEP,REJECT.
45 COMMON /REJECT/ AVNREJ,
46 * ALTMIN,ANEXP,THICKA,THICKD,CUTLN,EONCUT,
47 * FNPRIM
48 DOUBLE PRECISION AVNREJ(10)
49 REAL ALTMIN(10),ANEXP(10),THICKA(10),THICKD(10),
50 * CUTLN,EONCUT
51 LOGICAL FNPRIM
52*KEEP,RUNPAR.
53 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
54 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
55 * MONIOU,MDEBUG,NUCNUC,
56 * CETAPE,
57 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
58 * N1STTR,MDBASE,
59 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
60 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
61 * ,GHEISH,GHESIG
62 COMMON /RUNPAC/ DSN,HOST,USER
63 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
64 REAL STEPFC
65 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
66 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
67 * N1STTR,MDBASE
68 INTEGER CETAPE
69 CHARACTER*79 DSN
70 CHARACTER*20 HOST,USER
71
72 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
73 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
74 * ,GHEISH,GHESIG
75*KEEP,STACKE.
76 COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
77 DOUBLE PRECISION E(60),TIME(60)
78 REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
79 INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
80*KEND.
81 COMMON/THRESH/RMT2,RMSQ,ESCD2,AP,API,AE,UP,UE,TE,THMOLL
82 DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC
83 COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL
84 COMMON/ACLOCK/NCLOCK,JCLOCK
85 DOUBLE PRECISION THICK
86 CHARACTER MEDARR*24
87 DATA MEDARR/'AIR-NTP '/
88C---------------------------------------------------------------------
89
90 IF((DEBUG))WRITE(MDEBUG,* )'EGS4 :'
91C*** CONVERSION GEV --> MEV
92 E(1)= EEIN*1000.D0
93C*** CHECK ENERGY RANGE
94 IQ(1)=NINT(SECPAR(1))
95 IF ( IQ(1) .EQ. 1 ) THEN
96 IF ( E(1) .GT. UP ) THEN
97 CALL AUSGB2
98 WRITE(KMPO,91) EEIN
9991 FORMAT(' EGS4 : ENERGY OF GAMMA =',1P,E10.3,' GEV TOO HIGH')
100 STOP
101 ENDIF
102 ELSE
103 IF ( E(1) .GT. UE ) THEN
104 CALL AUSGB2
105 WRITE(KMPO,92) EEIN
10692 FORMAT(' EGS4 :ENERGY OF ELECTRON/POSITRON =',1P,E10.3,
107 * ' GEV TOO HIGH')
108 STOP
109 ENDIF
110 ENDIF
111C*** FILL IN STARTING COORDINATES
112 NP=1
113 TIME(1)=SECPAR(6)
114 X(1)=SECPAR(7)
115 Y(1)=-SECPAR(8)
116C*** STARTS IN HEIGHT 'Z' DOWNWARDS
117 Z(1)=-SECPAR(5)
118 IF (LLONGI) LPCTE(1)=MIN(NSTEP,INT(THICK(SECPAR(5))*THSTPI)+1)
119 SITHET=SQRT(1.D0-SECPAR(3)**2)
120C*** START DIRECTION COSINES
121 U(NP)=SITHET*COS(-SECPAR(4))
122 V(NP)=SITHET*SIN(-SECPAR(4))
123 W(NP)=SECPAR(3)
124C*** RENORMALIZATION OF DIRECTION COSINES
125 RADINV=1.5-0.5*(U(NP)**2+V(NP)**2+W(NP)**2)
126 U(NP)=U(NP)*RADINV
127 V(NP)=V(NP)*RADINV
128 W(NP)=W(NP)*RADINV
129 DNEAR(1)=0.
130 IGEN(1)=GEN
131 DO 101 K=1,5
132C *** DETERMINE START REGION
133 IF (-BOUND(K).LE.Z(1) .AND. -BOUND(K+1).GT.Z(1)) THEN
134 IR(1)=K+1
135 GO TO 110
136 END IF
137101 CONTINUE
138102 CONTINUE
139 CALL AUSGB2
140 WRITE(KMPO,120) -Z(1)*0.00001
141120 FORMAT (' EGS4 : START VALUE OF Z=',1P,E10.3,' KM NOT IN ATMOSPH
142 *ERE')
143 STOP
144110 CONTINUE
145 DO 111 IDET=1,NOBSLV
146C *** DETERMINE NEXT OBSERVATION LEVEL
147 IF (-Z(1).GE.OBSLVL(IDET)) THEN
148 IOBS(1)=IDET
149 GO TO 130
150 END IF
151111 CONTINUE
152112 CONTINUE
153 CALL AUSGB2
154 WRITE(KMPO,140) -Z(1)*0.01,OBSLVL(NOBSLV)*0.01
155140 FORMAT(' EGS4 : START VALUE OF Z= ',E10.3, ' M BELOW LOWEST DET'
156 *, 'ECTOR AT',E10.3,' M')
157 STOP
158C*** NEWOBS IS THE NEXT OBSERVATION LEVEL
159130 NEWOBS=IOBS(NP)
160 CALL SHOWER
161 IF(DEBUG) WRITE(MDEBUG,*) 'EGS4 : EGS-STACK EMPTY, EXIT'
162 RETURN
163 END
Note: See TracBrowser for help on using the repository browser.