source: branches/start/MagicSoft/Simulation/Corsika/Mmcs/compt.f@ 10113

Last change on this file since 10113 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: 4.2 KB
Line 
1 SUBROUTINE COMPT
2C VERSION 4.00 -- 26 JAN 1986/1900
3C******************************************************************
4C BUTCHER AND MESSEL'S CROSS SECTION EXPRESSION IS USED
5C (BUTCHER AND MESSEL, OP.CIT., P. 17-19,25), BUT THE
6C 1/EPSILON PART IS NOT SAMPLED IN THE WAY THAT THEY DO.
7C THIS ROUTINE CALLS THEIR 'EPSILON' VARIABLE BY THE NAME 'BR'.
8C BR=FINAL PHOTON ENERGY /INITIAL PHOTON ENERGY.
9C BR0 = MINIMUM BR = 1./(1.+2.*(E(NP)/RM))
10C MAXIMUM BR IS 1.
11C BUTCHER AND MESSEL'S EXPRESSION FOR THE DIFFERENTIAL CROSS
12C SECTION IS PROPORTIONAL TO
13C (1./BR+BR)*(1.-BR*SINTHE**2/(1.+BR*BR))
14C WE SHALL SAMPLE FROM THE FIRST FACTOR FROM THE INTERVAL (BR0,1)
15C AND USE THE SECOND FACTOR AS A REJECTION FUNCTION.
16C******************************************************************
17 DOUBLE PRECISION PEIG,PESG,PESE
18*KEEP,RANDPA.
19 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
20 DOUBLE PRECISION FAC,U1,U2
21 REAL RD(3000)
22 INTEGER ISEED(103,10),NSEQ
23 LOGICAL KNOR
24*KEEP,RUNPAR.
25 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
26 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
27 * MONIOU,MDEBUG,NUCNUC,
28 * CETAPE,
29 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
30 * N1STTR,MDBASE,
31 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
32 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
33 * ,GHEISH,GHESIG
34 COMMON /RUNPAC/ DSN,HOST,USER
35 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
36 REAL STEPFC
37 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
38 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
39 * N1STTR,MDBASE
40 INTEGER CETAPE
41 CHARACTER*79 DSN
42 CHARACTER*20 HOST,USER
43
44 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
45 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
46 * ,GHEISH,GHESIG
47*KEEP,STACKE.
48 COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
49 DOUBLE PRECISION E(60),TIME(60)
50 REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
51 INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
52*KEND.
53 COMMON/THRESH/RMT2,RMSQ,ESCD2,AP,API,AE,UP,UE,TE,THMOLL
54 COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
55 DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC
56 COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL
57 COMMON/ACLOCK/NCLOCK,JCLOCK
58C_____IF (NCLOCK.GT.JCLOCK) THEN
59C______WRITE(MDEBUG,* )' COMPT: NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
60C______CALL AUSGB2
61C_____END IF
62 PEIG=E(NP)
63 EIG=PEIG
64 EGP=EIG*RMI
65 BR0I=1.+2.*EGP
66 BR0=1./BR0I
67 ALPH1=LOG(BR0I)
68 ALPH2=EGP*(BR0I+1.)*BR0*BR0
69 SUMALP = ALPH1+ALPH2
70371 CONTINUE
71 CALL RMMAR(RNNO15,1,2)
72 IF (ALPH1.GE.SUMALP*RNNO15) THEN
73 CALL RMMAR(RNNO16,1,2)
74 BR=EXP(ALPH1*RNNO16)*BR0
75 ELSE
76 CALL RMMAR(RD,2,2)
77 BRP=RD(1)
78 RNNO18=RD(2)
79 IF (EGP.GE.(EGP+1.)*RNNO18) THEN
80 CALL RMMAR(RNNO19,1,2)
81 BRP=MAX(BRP,RNNO19)
82 END IF
83 BR=((BR0I-1.)*BRP+1.)*BR0
84 END IF
85 ESG=BR*EIG
86 A1MIBR = 1.-BR
87 TEMP=RM*A1MIBR/ESG
88 SINTHE=MAX(0.0,TEMP*(2.0-TEMP))
89 CALL RMMAR(RNNO20,1,2)
90 IF(((1.0-RNNO20)*(1.0+BR*BR).GE.BR*SINTHE))GO TO372
91 GO TO 371
92372 CONTINUE
93 SINTHE=SQRT(SINTHE)
94 COSTHE=1.0-TEMP
95 PESG=ESG
96 PESE=PEIG-PESG+PRM
97 ESE=PESE
98 CALL UPHI(2,1)
99 NP=NP+1
100 PSQ=ESE*ESE-RMSQ
101 IF (PSQ.LE.0.0) THEN
102 COSTHE=0.
103 SINTHE=-1.
104 ELSE
105 COSTHE=(ESE+ESG)*A1MIBR/SQRT(PSQ)
106 SINTHE=-SQRT(MAX(0.0,1.0-COSTHE*COSTHE))
107 END IF
108 CALL UPHI(3,2)
109 IF (ESE.LE.ESG) THEN
110 IQ(NP)=3
111 E(NP)=PESE
112 E(NP-1)=PESG
113 ELSE
114 IQ(NP)=1
115 IQ(NP-1)=3
116 E(NP)=PESG
117 E(NP-1)=PESE
118 T=U(NP)
119 U(NP)=U(NP-1)
120 U(NP-1)=T
121 T=V(NP)
122 V(NP)=V(NP-1)
123 V(NP-1)=T
124 T=W(NP)
125 W(NP)=W(NP-1)
126 W(NP-1)=T
127 END IF
128 RETURN
129 END
Note: See TracBrowser for help on using the repository browser.