source: trunk/MagicSoft/Simulation/Corsika/Mmcs/avage.f@ 20115

Last change on this file since 20115 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: 7.1 KB
Line 
1 SUBROUTINE AVAGE
2
3C-----------------------------------------------------------------------
4C AVE(ERAGE) AGE
5C
6C CALCULATES AVERAGE AGE AS A FUNCTION OF RADIUS
7C THIS SUBROUTINE IS CALLED FROM MAIN
8C-----------------------------------------------------------------------
9
10 IMPLICIT NONE
11*KEEP,BUFFS.
12 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
13 INTEGER MAXBUF,MAXLEN
14 PARAMETER (MAXBUF=39*7)
15 PARAMETER (MAXLEN=12)
16 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
17 * RUNE(MAXBUF),DATAB(MAXBUF)
18 INTEGER LH
19 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
20 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
21 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
22*KEEP,ELABCT.
23 COMMON /ELABCT/ ELCUT
24 DOUBLE PRECISION ELCUT(4)
25*KEEP,NKGI.
26 COMMON /NKGI/ SEL,SELLG,STH,ZEL,ZELLG,ZSL,DIST,
27 * DISX,DISY,DISXY,DISYX,DLAX,DLAY,DLAXY,DLAYX,
28 * OBSATI,RADNKG,RMOL,TLEV,TLEVCM,IALT
29 DOUBLE PRECISION SEL(10),SELLG(10),STH(10),ZEL(10),ZELLG(10),
30 * ZSL(10),DIST(10),
31 * DISX(-10:10),DISY(-10:10),
32 * DISXY(-10:10,2),DISYX(-10:10,2),
33 * DLAX (-10:10,2),DLAY (-10:10,2),
34 * DLAXY(-10:10,2),DLAYX(-10:10,2),
35 * OBSATI(2),RADNKG,RMOL(2),TLEV(10),TLEVCM(10)
36 INTEGER IALT(2)
37*KEEP,NKGS.
38 COMMON /NKGS/ CZX,CZY,CZXY,CZYX,SAH,SL,ZNE
39 DOUBLE PRECISION CZX(-10:10,2),CZY(-10:10,2),CZXY(-10:10,2),
40 * CZYX(-10:10,2),SAH(10),SL(10),ZNE(10)
41*KEEP,RUNPAR.
42 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
43 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
44 * MONIOU,MDEBUG,NUCNUC,
45 * CETAPE,
46 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
47 * N1STTR,MDBASE,
48 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
49 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
50 * ,GHEISH,GHESIG
51 COMMON /RUNPAC/ DSN,HOST,USER
52 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
53 REAL STEPFC
54 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
55 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
56 * N1STTR,MDBASE
57 INTEGER CETAPE
58 CHARACTER*79 DSN
59 CHARACTER*20 HOST,USER
60
61 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
62 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
63 * ,GHEISH,GHESIG
64*KEND.
65
66 DOUBLE PRECISION AJ,BJ,CJ,DF(10),SJ(10),SLLG,TH,ZF
67 INTEGER I,ID,IL,IOL,J,K,L
68C-----------------------------------------------------------------------
69
70 IF ( DEBUG ) WRITE(MDEBUG,*) 'AVAGE :'
71
72 IF ( FPRINT ) WRITE(MONIOU,1110) SHOWNO,ELCUT(3),ELCUT(4)
73 1110 FORMAT (/' ---------- NKG - OUTPUT OF SHOWER NO ',I10,
74 * ' --------------------------------'/
75 * ' ELECTRON/PHOTON THRESHOLD AT ',F10.5,' /',F10.5,' GEV')
76
77C LOOP OVER ALL DISTANCES WHERE ELECTRON NUMBER IS CALCULATED
78 DO 302 K = 1,2
79 IF ( OBSATI(K) .GE. 0.D0 ) THEN
80 DO 301 ID = -10,10
81 DLAX (ID,K) = DLAX (ID,K) + CZX (ID,K)
82 DLAY (ID,K) = DLAY (ID,K) + CZY (ID,K)
83 DLAXY(ID,K) = DLAXY(ID,K) + CZXY(ID,K)
84 DLAYX(ID,K) = DLAYX(ID,K) + CZYX(ID,K)
85 301 CONTINUE
86 ENDIF
87 302 CONTINUE
88
89C CALCULATE LONGITUDINAL SHOWER DEVELOPMENT
90 DO 311 IL = 1,IALT(1)
91 IF ( SL(IL) .GT. 0.D0 ) THEN
92 SEL(IL) = SEL(IL) + SL(IL)
93 SLLG = LOG10(SL(IL))
94 SELLG(IL) = SELLG(IL) + SLLG
95 ZEL(IL) = ZEL(IL) + SL(IL)**2
96 ZELLG(IL) = ZELLG(IL) + SLLG**2
97 ZF = ZNE(IL) / SL(IL)
98 CALL AGE( ZF,TH )
99C AGE PARAMETERS AVERAGED ON ALL SUBCASCADES AT THIS LEVEL
100 SAH(IL) = TH
101 STH(IL) = STH(IL) + TH
102 ZSL(IL) = ZSL(IL) + TH**2
103 ELSE
104 SAH(IL) = 0.D0
105 ENDIF
106 EVTE(175+IL) = SL (IL)
107 EVTE(185+IL) = SAH(IL)
108 EVTE(215+IL) = TLEV(IL)
109 EVTE(225+IL) = TLEVCM(IL)
110 311 CONTINUE
111
112C PRINT LONGITUDINAL SHOWER DEVELOPMENT
113 IF ( FPRINT ) WRITE(MONIOU,229)
114 * (I,TLEV(I),TLEVCM(I),SL(I),SAH(I),I=1,IALT(1))
115 229 FORMAT(
116 * /' LEVEL',2X,'THICKNESS',8X,'HEIGHT',5X,'ELECT. NUMBER',7X,'AGE'
117 * /' NO. ',2X,' G/CM**2',8X,' CM'/
118 * (' ',I4,F12.0,2X,F12.0,1X,F17.3,F10.3) )
119
120 DO 312 IOL = 1,2
121 IF ( OBSATI(IOL) .LT. 0.D0 ) GOTO 312
122C DETERMINE LOCAL AGE PARAMETER
123 DO 50 J = 1,9
124 IF ( CZX(J+1,IOL).GT.0.D0 .AND. CZX(-J-1,IOL).GT.0.D0 .AND.
125 * CZXY(J+1,IOL).GT.0.D0 .AND. CZXY(-J-1,IOL).GT.0.D0 .AND.
126 * CZYX(J+1,IOL).GT.0.D0 .AND. CZYX(-J-1,IOL).GT.0.D0 .AND.
127 * CZY(J+1,IOL).GT.0.D0 .AND. CZY(-J-1,IOL).GT.0.D0 ) THEN
128 AJ = 0.125D0 * (
129 * CZX(J,IOL) /CZX(J+1,IOL) + CZX(-J,IOL) /CZX(-J-1,IOL)
130 * + CZXY(J,IOL)/CZXY(J+1,IOL)+ CZXY(-J,IOL)/CZXY(-J-1,IOL)
131 * + CZYX(J,IOL)/CZYX(J+1,IOL)+ CZYX(-J,IOL)/CZYX(-J-1,IOL)
132 * + CZY(J,IOL) /CZY(J+1,IOL) + CZY(-J,IOL) /CZY(-J-1,IOL) )
133 ELSE
134 AJ = 0.D0
135 ENDIF
136 IF ( AJ .GT. 0.D0 ) THEN
137 BJ = DIST(J) / DIST(J+1)
138 CJ = (DIST(J)+RMOL(IOL)) / (DIST(J+1)+RMOL(IOL))
139 SJ(J) = LOG(AJ * BJ**2 * CJ**4.5D0) / LOG(BJ * CJ)
140 DF(J) = 0.5D0 * (DIST(J) + DIST(J+1))
141 ELSE
142 SJ(J) = 0.D0
143 DF(J) = 0.D0
144 ENDIF
145 50 CONTINUE
146
147 DO L = 1,10
148 EVTE(165+IOL*40+L) = SJ(L)
149 ENDDO
150
151 IF ( FPRINT ) THEN
152C WRITE LOCAL AGE PARAMETER
153 WRITE(MONIOU,60) IOL,OBSATI(IOL), (I,DF(I),SJ(I),I=1,9)
154 60 FORMAT(/' RADIAL BIN DISTANCE(CM) LOCAL AGE AT LEVEL NO.',
155 * I4,' AT HEIGHT:',F10.0,' CM'/
156 * (' ',I10,' ',F10.0,' ',F10.3 ) )
157
158C PRINT LATERAL ELECTRON DISTRIBUTION
159 WRITE(MONIOU,507) IOL,OBSATI(IOL)
160 507 FORMAT(/' LATERAL ELECTRON DENSITY (/CM**2) AT LEVEL NO.',
161 * I4,' AT HEIGHT:',F10.0,' CM'/
162 * ' --------------------------------------------------',
163 * '---------------------------'/
164 * ' DIST (CM) CZX CZY ',
165 * ' CZXY CZYX ')
166 WRITE(MONIOU,508) (DISX(I),CZX(I,IOL),CZY(I,IOL),
167 * CZXY(I,IOL),CZYX(I,IOL),I=-10,10)
168 508 FORMAT(' ',0P,F10.0,1P,4E15.5)
169 ENDIF
170
171 312 CONTINUE
172
173 DO L = 1,10
174 EVTE(195+L) = DIST(L)
175 EVTE(235+L) = DF(L)
176 ENDDO
177
178C WRITE NKG - SHOWER INFORMATION TO EVENT END BLOCK
179 DO 353 L = 1,21
180 EVTE( 7+L) = CZX (-11+L,1)
181 EVTE( 28+L) = CZY (-11+L,1)
182 EVTE( 49+L) = CZXY(-11+L,1)
183 EVTE( 70+L) = CZYX(-11+L,1)
184 EVTE( 91+L) = CZX (-11+L,2)
185 EVTE(112+L) = CZY (-11+L,2)
186 EVTE(133+L) = CZXY(-11+L,2)
187 EVTE(154+L) = CZYX(-11+L,2)
188 353 CONTINUE
189
190 RETURN
191 END
Note: See TracBrowser for help on using the repository browser.