source: trunk/MagicSoft/Simulation/Corsika/Mmcs/photon.f@ 402

Last change on this file since 402 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: 10.8 KB
Line 
1 SUBROUTINE PHOTON(IRCODE)
2C VERSION 4.00 -- 26 JAN 1986/1900
3C******************************************************************
4 DOUBLE PRECISION PEIG
5 COMMON/BOUNDS/ECUT(6),PCUT(6),VACDST
6*KEEP,BUFFS.
7 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
8 INTEGER MAXBUF,MAXLEN
9 PARAMETER (MAXBUF=39*7)
10 PARAMETER (MAXLEN=12)
11 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
12 * RUNE(MAXBUF),DATAB(MAXBUF)
13 INTEGER LH
14 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
15 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
16 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
17*KEEP,CEREN1.
18 COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
19 * CERSIZ,LCERFI
20 DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
21 REAL CERSIZ
22 LOGICAL LCERFI
23*KEEP,EPCONT.
24 COMMON/EPCONT/ EDEP,RATIO,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,IDISC,
25 * IROLD,IRNEW,RHOFAC, EOLD,ENEW,EKE,ELKE,BETA2,GLE,
26 * TSCAT,IAUSFL
27 DOUBLE PRECISION EDEP,RATIO
28 REAL TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,RHOFAC,EOLD,ENEW,
29 * EKE,ELKE,BETA2,GLE,TSCAT
30 INTEGER IDISC,IROLD,IRNEW,IAUSFL(29)
31*KEND.
32 COMMON/GEOM/ZALTIT,BOUND(6),NEWOBS,OBSLVL(10)
33*KEEP,LONGI.
34 COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
35 * NSTEP,LLONGI,FLGFIT
36 DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
37 * SPLONG(0:1040,9),THSTEP,THSTPI
38 INTEGER NSTEP
39 LOGICAL LLONGI,FLGFIT
40*KEND.
41 COMMON /MEDIA/ NMED, RLC,RLDU,RLDUI,RHO,MSGE,MGE,MSEKE,MEKE,MLEKE,
42 *MCMFP,MRANGE,IRAYLM,HBARO(6),HBAROI(6)
43 CHARACTER MEDIA*24
44 COMMON/MEDIAC/MEDIA
45 COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(6),RHOR(6),IRAYLR(6)
46 DOUBLE PRECISION PRRMMU
47 COMMON/MUON/PRRMMU,RMMU,RMMUT2
48*KEEP,OBSPAR.
49 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
50 * THETPR,PHIPR,NOBSLV
51 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
52 * THETAP,THETPR(2),PHIP,PHIPR(2)
53 INTEGER NOBSLV
54*KEEP,PARPAR.
55 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
56 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
57 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
58 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
59 INTEGER ITYPE,LEVL
60*KEND.
61 COMMON/PHOTIN/EBINDA,GE0,GE1, MPGEM(1),GMFP0(500),GMFP1(500),GBR10
62 *(500),GBR11(500),GBR20(500),GBR21(500),GBR30(500),GBR31(500),GBR40
63 *(500),GBR41(500),NGR,RCO0,RCO1, RSCT0(100),RSCT1(100), COHE0(500),
64 *COHE1(500)
65 DOUBLE PRECISION PI0MSQ
66 COMMON/PION/PI0MSQ,PITHR,PICMAS,PI0MAS,AMASK0,AMASKC,AMASPR,AMASNT
67 *
68*KEEP,RANDPA.
69 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
70 DOUBLE PRECISION FAC,U1,U2
71 REAL RD(3000)
72 INTEGER ISEED(103,10),NSEQ
73 LOGICAL KNOR
74*KEEP,REJECT.
75 COMMON /REJECT/ AVNREJ,
76 * ALTMIN,ANEXP,THICKA,THICKD,CUTLN,EONCUT,
77 * FNPRIM
78 DOUBLE PRECISION AVNREJ(10)
79 REAL ALTMIN(10),ANEXP(10),THICKA(10),THICKD(10),
80 * CUTLN,EONCUT
81 LOGICAL FNPRIM
82*KEEP,RUNPAR.
83 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
84 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
85 * MONIOU,MDEBUG,NUCNUC,
86 * CETAPE,
87 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
88 * N1STTR,MDBASE,
89 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
90 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
91 * ,GHEISH,GHESIG
92 COMMON /RUNPAC/ DSN,HOST,USER
93 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
94 REAL STEPFC
95 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
96 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
97 * N1STTR,MDBASE
98 INTEGER CETAPE
99 CHARACTER*79 DSN
100 CHARACTER*20 HOST,USER
101
102 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
103 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
104 * ,GHEISH,GHESIG
105*KEEP,STACKE.
106 COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
107 DOUBLE PRECISION E(60),TIME(60)
108 REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
109 INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
110*KEND.
111 COMMON/THRESH/RMT2,RMSQ,ESCD2,AP,API,AE,UP,UE,TE,THMOLL
112 COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
113 DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC
114 COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL
115 COMMON/ACLOCK/NCLOCK,JCLOCK
116 DOUBLE PRECISION THICK
117
118
119
120C_____IF (NCLOCK.GT.JCLOCK) THEN
121C______WRITE(MDEBUG,* )' PHOTON:NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
122C______CALL AUSGB2
123C_____END IF
124 NEWOBS=IOBS(NP)
125 IRCODE=1
126 PEIG=E(NP)
127 EIG=PEIG
128 IRL=IR(NP)
129 MEDIUM=MED(IRL)
130 IF((EIG.LE.PCUT(IRL)))GO TO 970
131980 CONTINUE
132981 CONTINUE
133 GLE=LOG(EIG)
134 CALL RMMAR(RNNO35,1,2)
135 IF ((RNNO35.EQ.0.0)) THEN
136 RNNO35=1.E-30
137 END IF
138 DPMFP=-ALOG(RNNO35)
139 IROLD=IR(NP)
1401030 CONTINUE
1411031 CONTINUE
142 IF (MEDIUM.NE.0) THEN
143 LGLE=GE1*GLE+GE0
144 GMFPR0=GMFP1(LGLE)*GLE+GMFP0(LGLE)
145 END IF
1461040 CONTINUE
1471041 CONTINUE
148 IF (MEDIUM.EQ.0) THEN
149 TSTEP=VACDST
150 ELSE
151 RHOFAC=RHOR(IRL)/RHO
152 RHOFI=1./RHOFAC
153 GMFP=GMFPR0*RHOFI
154 IF ((IRAYLR(IRL).EQ.1)) THEN
155 COHFAC=COHE1(LGLE)*GLE+COHE0(LGLE)
156 GMFP=GMFP*COHFAC
157 END IF
158 TSTEP=GMFP*DPMFP
159 ALTEXP=EXP(-Z(NP)*HBAROI(IRL))
160 TSTEP=TSTEP*ALTEXP
161 DISC=W(NP)*TSTEP*HBAROI(IRL)
162 IF (ABS(DISC).LT.0.065) THEN
163 TSTEP=TSTEP*(1.-0.5*DISC*(1.-0.6666667*DISC* (1.-0.75*DISC *
164 * (1.-0.8*DISC))))
165 ELSE IF(DISC.LE.-1.) THEN
166 TSTEP=VACDST
167 ELSE
168 TSTEP=TSTEP/DISC*LOG(DISC+1.)
169 END IF
170 END IF
171 IRNEW=IR(NP)
172 IDISC=0
173 USTEP=TSTEP
174 TUSTEP=USTEP
175 IF((USTEP.GT.DNEAR(NP)))CALL HOWFAR
176 IF((IDISC.GT.0))GO TO 1000
177 VSTEP=USTEP
178 TVSTEP=VSTEP
179 EDEP=PZERO
180 USTEPU=USTEP
181 DISC=W(NP)*USTEPU*HBAROI(IRL)
182 USTEPU=USTEPU/ALTEXP
183 IF (ABS(DISC).LT.0.16) THEN
184 USTEPU=USTEPU*(1.+.5*DISC*(1.+.33333333*DISC* (1.+0.25*DISC* (
185 * 1.+0.2*DISC))))
186 ELSE
187 USTEPU=USTEPU/DISC*(EXP(DISC)-1.)
188 END IF
189 X(NP)=X(NP)+U(NP)*USTEP
190 Y(NP)=Y(NP)+V(NP)*USTEP
191 ZOLD =Z(NP)
192 Z(NP)=Z(NP)+W(NP)*USTEP
193 TIME(NP)=TIME(NP)+TVSTEP*VC
194
195
196C ADD PHOTONS TO THE LONGITUDINAL DEVELOPMENT
197 IF ( LLONGI ) THEN
198C FIND FIRST THE EQUIVALENT LEVELS
199C IF STARTING POINT BELOW LOWEST LEVEL THEN DON'T CHECK
200 IF ( HLONG(NSTEP) .LE. -ZOLD ) THEN
201 LPCT1 = LPCTE(NP)
202C Z NEW IS PROBABLY ONLY LITTLE BELOW Z OLD, THEREFORE INCREMENTAL SEARCH
203 DO 6002 I1 = LPCT1,NSTEP
204 IF ( HLONG(I1) .LT. -Z(NP) ) GOTO 6003
205 6002 CONTINUE
206 I1 = NSTEP + 1
207 6003 CONTINUE
208 LPCT2 = I1 - 1
209 DO 485 I=LPCT1,LPCT2
210 PLONG(I,1) = PLONG(I,1) + 1.D0
211 485 CONTINUE
212 LPCTE(NP) = LPCT2 + 1
213 ENDIF
214 ENDIF
215 DNEAR(NP)=DNEAR(NP)-USTEP
216 IF (MEDIUM.NE.0) THEN
217 DPMFP=MAX(0.,DPMFP-USTEPU/GMFP)
218 END IF
219 IROLD=IR(NP)
220 MEDOLD=MEDIUM
221 IF (IRNEW.NE.IROLD) THEN
222 IR(NP)=IRNEW
223 IRL=IRNEW
224 MEDIUM=MED(IRL)
225 IF((EIG.LE.PCUT(IRL)))GO TO 970
226 END IF
227 IF (NEWOBS.GT.IOBS(NP)) THEN
228 CALL AUSGAB
229 IOBS(NP)=NEWOBS
230 END IF
231 IF((IDISC.LT.0))GO TO 1000
232 IF((MEDIUM.NE.MEDOLD))GO TO 1042
233 IF((MEDIUM.NE.0.AND.DPMFP.LE.1.E-6))GO TO 1032
234 GO TO 1041
2351042 CONTINUE
236 GO TO 1031
2371032 CONTINUE
238 IF ((IRAYLR(IRL).EQ.1)) THEN
239 CALL RMMAR(RNNO37,1,2)
240 IF ((RNNO37.LE.(1.0-COHFAC))) THEN
2411050 CONTINUE
2421051 CONTINUE
243 CALL RMMAR(XXX,1,2)
244 LXXX=RCO1*XXX+RCO0
245 X2=RSCT1(LXXX)*XXX+RSCT0(LXXX)
246 Q2=X2*RMSQ*.23547885E-02
247 COSTHE=1.-Q2/(2.*E(NP)*E(NP))
248 IF((ABS(COSTHE).GT.1.0))GO TO 1050
249 CSQTHE=COSTHE*COSTHE
250 REJF=(1.0+CSQTHE)*.5
251 CALL RMMAR(RNNORJ,1,2)
252 IF((RNNORJ.LE.REJF))GO TO1052
253 GO TO 1051
2541052 CONTINUE
255 SINTHE=SQRT(AMAX1(0.,1.0-CSQTHE))
256 CALL UPHI(2,1)
257 GOTO 980
258 END IF
259 END IF
260 IF ( .NOT. FNPRIM ) THEN
261 X(1)=0.
262 Y(1)=0.
263 EVTH(5)=X(1)
264 EVTH(6)=-Y(1)
265 IF (FIX1I) THEN
266 Z(1)=-FIXHEI
267 NP=1
268 LPCTE(1)=MIN(NSTEP,INT(THICK(FIXHEI)*THSTPI)+1)
269 SITHET=SQRT(1.D0-SECPAR(3)**2)
270 U(1)=SITHET*COS(-SECPAR(4))
271 V(1)=SITHET*SIN(-SECPAR(4))
272 W(1)=SECPAR(3)
273 RADINV=1.5-0.5*(U(1)**2+V(1)**2+W(1)**2)
274 U(1)=U(1)*RADINV
275 V(1)=V(1)*RADINV
276 W(1)=W(1)*RADINV
277 END IF
278 EVTH(7)=-Z(1)
279 CALL TOBUF(EVTH,0)
280 IF (LCERFI) CALL TOBUFC(EVTH,0)
281 CALL COORIN(DBLE(-Z(1)))
282 TIME(1)=0.D0
283 FNPRIM =.TRUE.
284 IF (FPRINT) THEN
285 WRITE(KMPO,* )' FIRST INTERACTION AT ',EVTH(7)*0.01,' M'
286 END IF
287 END IF
288 CALL RMMAR(RNNO36,1,2)
289 GBR1=GBR11(LGLE)*GLE+GBR10(LGLE)
290 IF ((RNNO36.LE.GBR1).AND.(E(NP).GT.RMT2)) THEN
291 CALL PAIR
292 GO TO 982
293 END IF
294 GBR2=GBR21(LGLE)*GLE+GBR20(LGLE)
295 IF (RNNO36.LT.GBR2) THEN
296 CALL COMPT
297 IF((IQ(NP).NE.1))GO TO 982
298 GO TO1060
299 END IF
300 GBR4=GBR41(LGLE)*GLE+GBR40(LGLE)
301 IF (RNNO36.GE.GBR4 .AND. E(NP).GT.RMMUT2) THEN
302 CALL MUPAIR
303 GO TO 982
304 END IF
305 GBR3=GBR31(LGLE)*GLE+GBR30(LGLE)
306 IF (RNNO36.GE.GBR3 .AND. E(NP).GT.PITHR) THEN
307 CALL PIGEN
308 IF (NP.EQ.0) THEN
309 IRCODE=2
310 RETURN
311 END IF
312 GO TO 982
313 ELSE
314 CALL PHOTO
315 IF (NP.EQ.0) THEN
316 IRCODE=2
317 RETURN
318 END IF
319 IF((IQ(NP).EQ.3))GO TO 982
320 END IF
3211060 PEIG=E(NP)
322 EIG=PEIG
323 IF((EIG.LT.PCUT(IRL)))GO TO 970
324 GO TO 981
325982 CONTINUE
326 RETURN
327970 IF (EIG.GT.AP) THEN
328 IDR=1
329 ELSE
330 IDR=2
331 END IF
332 EDEP=PEIG
333 IRCODE=2
334 NP=NP-1
335 RETURN
3361000 EDEP=PEIG
337 IRCODE=2
338 NP=NP-1
339 RETURN
340 END
Note: See TracBrowser for help on using the repository browser.