source: branches/start/MagicSoft/Simulation/Corsika/Mmcs/hatch.f@ 18569

Last change on this file since 18569 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: 13.2 KB
Line 
1 SUBROUTINE HATCH
2C VERSION 4.00 -- 26 JAN 1986/1900
3C******************************************************************
4C SETUP WHICH THE USER IS EXPECTED TO DO BEFORE CALLING HATCH IS:
5C 1. SET 'NMED' TO THE NUMBER OF MEDIA TO BE USED.
6C 2. SET THE ARRAY 'MEDIA', WHICH CONTAINS THE NAMES OF THE
7C MEDIA THAT ARE DESIRED. THE CHARACTER FORMAT IS A1, SO
8C THAT MEDIA(IB,IM) CONTAINS THE IB'TH BYTE OF THE NAME OF
9C THE IM'TH MEDIUM IN A1 FORMAT.
10C 3. SET 'DUNIT', THE DISTANCE UNIT TO BE USED.
11C DUNIT.GT.0 MEANS VALUE OF DUNIT IS LENGTH OF DISTANCE UNIT
12C CENTIMETERS. DUNIT.LT.0 MEANS USE THE RADIATION LENGTH OF
13C THE ABS(DUNIT)'TH MEDIUM FOR THE DISTANCE UNIT.
14C 4. FILL THE ARRAY 'MED' WITH THE MEDIUM INDICES FOR THE
15C REGIONS.
16C 5. FILL ARRAYS 'ECUT' AND 'PCUT' WITH THE ELECTRON AND PHOTON
17C CUT-OFF ENERGIES FOR EACH REGION RESPECTIVELY. SETUP WILL
18C RAISE THESE IF NECESSARY TO MAKE THEM AT LEAST AS LARGE AS
19C THE REGION'S MEDIUM'S AE AND AP RESPECTIVELY.
20C 6. FILL 'MED' ARRAY. MED(IR) IS THE MEDIUM INDEX FOR REGION
21C IR. A ZERO MEDIUM INDEX MEANS THE REGION IS IN A VACUUM.
22C 7. FILL THE ARRAY 'IRAYLR' WITH 1 FOR EACH REGION IN WHICH
23C RAYLEIGH (COHERENT) SCATTERING IS TO BE INCLUDED.
24C******************************************************************
25 CHARACTER MBUF*72,MDLABL*8
26 DIMENSION ZEROS(3)
27CNOTE: ABOVE IS ZEROS OF SINE, 0,PI,TWOPI
28 COMMON/BOUNDS/ECUT(6),PCUT(6),VACDST
29 COMMON/BREMPR/DL1(6),DL2(6),DL3(6),DL4(6),DL5(6),DL6(6),DELCM, ALP
30 *HI(2),BPAR(2),DELPOS(2),PWR2I(50)
31 COMMON/ELECIN/EKELIM,ICOMP,EKE0,EKE1,CMFP0,CMFP1,RANGE0,RANGE1, XR
32 *0,TEFF0,BLCC,XCC,PICMP0(1),PICMP1(1),EICMP0(1),EICMP1(1),MPEEM(1),
33 * ESIG0(500),ESIG1(500),PSIG0(500),PSIG1(500),EDEDX0(500),EDEDX1(50
34 *0),PDEDX0(500),PDEDX1(500),EBR10(500),EBR11(500),PBR10(500),PBR11(
35 *500),PBR20(500),PBR21(500),TMXS0(500),TMXS1(500),CMFPE0(1),CMFPE1(
36 *1),CMFPP0(1),CMFPP1(1),ERANG0(1),ERANG1(1),PRANG0(1),PRANG1(1),CXC
37 *2E0(1),CXC2E1(1),CXC2P0(1),CXC2P1(1),CLXAE0(1),CLXAE1(1),CLXAP0(1)
38 *,CLXAP1(1), THR0(1,1),THR1(1,1),THR2(1,1),THRI0(1,1),THRI1(1,1),TH
39 *RI2(1,1),FSTEP(16),FSQR(16),MSMAP(200), VERT1(1000),VERT2(100,16),
40 *MSTEPS,JRMAX,MXV1, MXV2,NBLC,NRNTH,NRNTHI,BLC0,BLC1,RTHR0,RTHR1,RT
41 *HRI0,RTHRI1
42 COMMON /MEDIA/ NMED, RLC,RLDU,RLDUI,RHO,MSGE,MGE,MSEKE,MEKE,MLEKE,
43 *MCMFP,MRANGE,IRAYLM,HBARO(6),HBAROI(6)
44 CHARACTER MEDIA*24
45 COMMON/MEDIAC/MEDIA
46 COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(6),RHOR(6),IRAYLR(6)
47 COMMON/PHOTIN/EBINDA,GE0,GE1, MPGEM(1),GMFP0(500),GMFP1(500),GBR10
48 *(500),GBR11(500),GBR20(500),GBR21(500),GBR30(500),GBR31(500),GBR40
49 *(500),GBR41(500),NGR,RCO0,RCO1, RSCT0(100),RSCT1(100), COHE0(500),
50 *COHE1(500)
51*KEEP,RANDPA.
52 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
53 DOUBLE PRECISION FAC,U1,U2
54 REAL RD(3000)
55 INTEGER ISEED(103,10),NSEQ
56 LOGICAL KNOR
57*KEEP,RUNPAR.
58 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
59 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
60 * MONIOU,MDEBUG,NUCNUC,
61 * CETAPE,
62 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
63 * N1STTR,MDBASE,
64 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
65 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
66 * ,GHEISH,GHESIG
67 COMMON /RUNPAC/ DSN,HOST,USER
68 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
69 REAL STEPFC
70 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
71 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
72 * N1STTR,MDBASE
73 INTEGER CETAPE
74 CHARACTER*79 DSN
75 CHARACTER*20 HOST,USER
76
77 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
78 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
79 * ,GHEISH,GHESIG
80*KEEP,STACKE.
81 COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
82 DOUBLE PRECISION E(60),TIME(60)
83 REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
84 INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
85*KEND.
86 COMMON/THRESH/RMT2,RMSQ,ESCD2,AP,API,AE,UP,UE,TE,THMOLL
87 COMMON/UPHIIN/SINC0,SINC1,SIN0(20002),SIN1(20002)
88 COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
89 DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC
90 COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL
91 DATA MDLABL/' MEDIUM='/,LMDL/8/,LMDN/24/,DUNITO/1./
92 DATA I1ST/1/,NSINSS/37/,MXSINC/20002/,ISTEST/0/,NRNA/1000/
93510 FORMAT(1X,14I5)
94520 FORMAT(1X,1P,5E14.5)
95530 FORMAT(A72)
96 IF (I1ST.NE.0) THEN
97 I1ST=0
98 PRM=RM
99 RMI=1./PRM
100 PRMT2=2.D0*PRM
101 PZERO=0.0D0
102 NISUB=MXSINC-2
103 FNSSS=NSINSS
104 WID=PI5D2/REAL(NISUB)
105 WSS=WID/(FNSSS-1.0)
106 ZEROS(1)=0.
107 ZEROS(2)=PI
108 ZEROS(3)=TWOPI
109 DO 541 ISUB=1,MXSINC
110 SX=0.
111 SY=0.
112 SXX=0.
113 SXY=0.
114 XS0=WID*REAL(ISUB-2)
115 XS1=XS0+WID
116 IZ=0
117 DO 551 IZZ=1,3
118 IF ((XS0.LE.ZEROS(IZZ)).AND.(ZEROS(IZZ).LE.XS1)) THEN
119 IZ=IZZ
120 GO TO552
121 END IF
122551 CONTINUE
123552 CONTINUE
124 IF (IZ.EQ.0) THEN
125 XSI=XS0
126 ELSE
127 XSI=ZEROS(IZ)
128 END IF
129 DO 561 ISS=1,NSINSS
130 XS=WID*REAL(ISUB-2)+WSS*REAL(ISS-1)-XSI
131 YS=SIN(XS+XSI)
132 SX=SX+XS
133 SY=SY+YS
134 SXX=SXX+XS*XS
135 SXY=SXY+XS*YS
136561 CONTINUE
137562 CONTINUE
138 IF (IZ.NE.0) THEN
139 SIN1(ISUB)=SXY/SXX
140 SIN0(ISUB)=-SIN1(ISUB)*XSI
141 ELSE
142 DEL=FNSSS*SXX-SX*SX
143 SIN1(ISUB)=(FNSSS*SXY-SY*SX)/DEL
144 SIN0(ISUB)=(SY*SXX-SX*SXY)/DEL - SIN1(ISUB)*XSI
145 END IF
146541 CONTINUE
147542 CONTINUE
148 SINC0=2.0
149 SINC1=1.0/WID
150 IF (ISTEST.NE.0) THEN
151 ADEV=0.
152 RDEV=0.
153 S2C2MN=10.
154 S2C2MX=0.
155 DO 571 ISUB=1,NISUB
156 DO 581 ISS=1,NSINSS
157 THETA=WID*REAL(ISUB-1)+WSS*REAL(ISS-1)
158 CTHET=PI5D2-THETA
159 LTHETA=SINC1*THETA+SINC0
160 LCTHET=SINC1*CTHET+SINC0
161 SINTHE=SIN1(LTHETA)*THETA+SIN0(LTHETA)
162 COSTHE=SIN1(LCTHET)*CTHET+SIN0(LCTHET)
163 SINT=SIN(THETA)
164 COST=COS(THETA)
165 ASD=ABS(SINTHE-SINT)
166 ACD=ABS(COSTHE-COST)
167 ADEV=MAX(ADEV,ASD,ACD)
168 IF((SINT.NE.0.0))RDEV=MAX(RDEV,ASD/ABS(SINT))
169 IF((COST.NE.0.0))RDEV=MAX(RDEV,ACD/ABS(COST))
170 S2C2=SINTHE**2+COSTHE**2
171 S2C2MN=MIN(S2C2MN,S2C2)
172 S2C2MX=MAX(S2C2MX,S2C2)
173 IF (ISUB.LT.11) THEN
174 WRITE(KMPO,590)THETA,SINTHE,SINT,COSTHE,COST
175590 FORMAT(1P,5E20.7)
176 END IF
177581 CONTINUE
178582 CONTINUE
179571 CONTINUE
180572 CONTINUE
181 WRITE(KMPO,600)MXSINC,NSINSS
182600 FORMAT(' SINE TESTS,MXSINC,NSINSS=',2I5)
183 WRITE(KMPO,610)ADEV,RDEV,S2C2MN,S2C2MX
184610 FORMAT(' ADEV,RDEV,S2C2(MN,MX) =',1P,4E16.8)
185 ADEV=0.
186 RDEV=0.
187 S2C2MN=10.
188 S2C2MX=0.
189 DO 621 IRN=1,NRNA
190 CALL RMMAR(THETA,1,2)
191 THETA=THETA*PI5D2
192 CTHET=PI5D2-THETA
193 LTHETA=SINC1*THETA+SINC0
194 LCTHET=SINC1*CTHET+SINC0
195 SINTHE=SIN1(LTHETA)*THETA+SIN0(LTHETA)
196 COSTHE=SIN1(LCTHET)*CTHET+SIN0(LCTHET)
197 SINT=SIN(THETA)
198 COST=COS(THETA)
199 ASD=ABS(SINTHE-SINT)
200 ACD=ABS(COSTHE-COST)
201 ADEV=MAX(ADEV,ASD,ACD)
202 IF((SINT.NE.0.0))RDEV=MAX(RDEV,ASD/ABS(SINT))
203 IF((COST.NE.0.0))RDEV=MAX(RDEV,ACD/ABS(COST))
204 S2C2=SINTHE**2+COSTHE**2
205 S2C2MN=MIN(S2C2MN,S2C2)
206 S2C2MX=MAX(S2C2MX,S2C2)
207621 CONTINUE
208622 CONTINUE
209 WRITE(KMPO,630)NRNA
210630 FORMAT(' TEST AT ',I7,' RANDOM ANGLES IN (0,5*PI/2)')
211 WRITE(KMPO,640)ADEV,RDEV,S2C2MN,S2C2MX
212640 FORMAT(' ADEV,RDEV,S2C2(MN,MX) =',1P,4E16.8)
213 END IF
214 P=1.
215 DO 651 I=1,50
216 PWR2I(I)=P
217 P=P*.5
218651 CONTINUE
219652 CONTINUE
220 END IF
221 DO 661 IM=1,NMED
222670 CONTINUE
223 DO 671 I=1,6
224 IF (IRAYLR(I).EQ.1.AND.MED(I).EQ.IM) THEN
225 IRAYLM=1
226 GO TO 672
227 END IF
228671 CONTINUE
229672 CONTINUE
230661 CONTINUE
231662 CONTINUE
232 REWIND KMPI
233 NM=0
234 DO 681 IM=1,NMED
235 LOK=0
236 IF (IRAYLM.EQ.1) THEN
237 WRITE(KMPO,690)IM
238690 FORMAT(' RAYLEIGH OPTION REQUESTED FOR MEDIUM NUMBER',I3,/)
239 END IF
240681 CONTINUE
241682 CONTINUE
242700 CONTINUE
243701 CONTINUE
244710 CONTINUE
245711 CONTINUE
246 READ(KMPI,530,END=720)MBUF
247 DO 731 IB=1,LMDL
248 IF((MBUF(IB:IB).NE.MDLABL(IB:IB)))GO TO 711
249731 CONTINUE
250732 CONTINUE
251740 CONTINUE
252 DO 741 IM=1,NMED
253 DO 751 IB=1,LMDN
254 IL=LMDL+IB
255 IF((MBUF(IL:IL).NE.MEDIA(IB:IB)))GO TO 741
256 IF((IB.EQ.LMDN))GO TO 712
257751 CONTINUE
258752 CONTINUE
259741 CONTINUE
260742 CONTINUE
261 GO TO 711
262712 CONTINUE
263 IF((LOK.NE.0))GO TO 710
264 LOK=1
265 NM=NM+1
266 WRITE(KMPO,760)IM,MBUF
267760 FORMAT(' DATA FOR MEDIUM #',I3,', WHICH IS:',A72)
268 READ(KMPI,770)(MBUF(I:I),I=1,5),RHO,NE
269770 FORMAT(5A1,5X,F11.0,4X,I2)
270 WRITE(KMPO,780)(MBUF(I:I),I=1,5),RHO,NE
271780 FORMAT(5A1,',RHO=',1P,G11.4, ',NE=',I2,',COMPOSITION IS :')
272 DO 791 IE=1,NE
273 READ(KMPI,530)MBUF
274 WRITE(KMPO,530)MBUF
275791 CONTINUE
276792 CONTINUE
277 READ(KMPI,520)RLC,AE,AP,UE,UP
278 TE=AE-RM
279 THMOLL=TE*2. + RM
280 READ(KMPI,510)MSGE,MGE,MSEKE,MEKE,MLEKE,MCMFP,MRANGE,IRAYL
281 NSGE=MSGE
282 NGE=MGE
283 NSEKE=MSEKE
284 NEKE=MEKE
285 NLEKE=MLEKE
286 NCMFP=MCMFP
287 NRANGE=MRANGE
288 READ(KMPI,520)(DL1(I),DL2(I),DL3(I),DL4(I),DL5(I),DL6(I),I=1,6)
289 READ(KMPI,520)DELCM,(ALPHI(I),BPAR(I),DELPOS(I),I=1,2)
290 READ(KMPI,520)XR0,TEFF0,BLCC,XCC
291 READ(KMPI,520)EKE0,EKE1
292 READ(KMPI,520)(ESIG0(I),ESIG1(I),PSIG0(I),PSIG1(I),EDEDX0(I),EDED
293 * X1(I),PDEDX0(I),PDEDX1(I),EBR10(I),EBR11(I),PBR10(I),PBR11(I),PBR
294 * 20(I),PBR21(I),TMXS0(I),TMXS1(I),I=1,NEKE)
295 READ(KMPI,520)EBINDA,GE0,GE1
296 READ(KMPI,520)(GMFP0(I),GMFP1(I),GBR10(I),GBR11(I),GBR20(I),GBR21
297 * (I),GBR30(I),GBR31(I),GBR40(I),GBR41(I),I=1,NGE)
298 IF (IRAYLM.EQ.1.AND.IRAYL.NE.1) THEN
299 WRITE(KMPO,800)IM
300800 FORMAT(' STOPPED IN HATCH: REQUESTED RAYLEIGH OPTION FOR MEDIUM'
301 * ,I3/ ' BUT RAYLEIGH DATA NOT INCLUDED IN DATA CREATED BY PEGS.')
302 STOP
303 END IF
304 IF (IRAYL.EQ.1) THEN
305 READ(KMPI,510)NGR
306 NGRIM=NGR
307 READ(KMPI,520)RCO0,RCO1
308 READ(KMPI,520)(RSCT0(I),RSCT1(I),I=1,NGRIM)
309 READ(KMPI,520)(COHE0(I),COHE1(I),I=1,NGE)
310 IF (IRAYLM.NE.1) THEN
311 WRITE(KMPO,810)IM
312810 FORMAT(' RAYLEIGH DATA AVAILABLE FOR MEDIUM',I3,' BUT OPTION ',
313 * 'NOT REQUESTED.',/)
314 END IF
315 END IF
316 IF((NM.GE.NMED))GO TO702
317 GO TO 701
318702 CONTINUE
319 DUNITR=DUNIT
320 IF (DUNIT.LT.0.0) THEN
321 MD=MAX(1,MIN(1,IFIX(-DUNIT)))
322 DUNIT=RLC
323 END IF
324 IF (DUNIT.NE.1.0) THEN
325 WRITE(KMPO,820)DUNITR,DUNIT
326820 FORMAT(' DUNIT REQUESTED&USED ARE:',1P,2E14.5,'(CM.)')
327 END IF
328 DO 831 IM=1,NMED
329 DFACT=RLC/DUNIT
330 DFACTI=1.0/DFACT
331 I=1
332 GO TO 843
333841 I=I+1
334843 IF(I-(MEKE).GT.0)GO TO 842
335 ESIG0(I)=ESIG0(I)*DFACTI
336 ESIG1(I)=ESIG1(I)*DFACTI
337 PSIG0(I)=PSIG0(I)*DFACTI
338 PSIG1(I)=PSIG1(I)*DFACTI
339 EDEDX0(I)=EDEDX0(I)*DFACTI
340 EDEDX1(I)=EDEDX1(I)*DFACTI
341 PDEDX0(I)=PDEDX0(I)*DFACTI
342 PDEDX1(I)=PDEDX1(I)*DFACTI
343 TMXS0(I)=TMXS0(I)*DFACT
344 TMXS1(I)=TMXS1(I)*DFACT
345 GO TO 841
346842 CONTINUE
347 I=1
348 GO TO 853
349851 I=I+1
350853 IF(I-(MLEKE).GT.0)GO TO 852
351 ERANG0(I)=ERANG0(I)*DFACT
352 ERANG1(I)=ERANG1(I)*DFACT
353 PRANG0(I)=PRANG0(I)*DFACT
354 PRANG1(I)=PRANG1(I)*DFACT
355 GO TO 851
356852 CONTINUE
357 TEFF0=TEFF0*DFACT
358 BLCC=BLCC*DFACTI
359 XCC=XCC*SQRT(DFACTI)
360 RLDU=RLC/DUNIT
361 RLDUI=1./RLDU
362 I=1
363 GO TO 863
364861 I=I+1
365863 IF(I-(MGE).GT.0)GO TO 862
366 GMFP0(I)=GMFP0(I)*DFACT
367 GMFP1(I)=GMFP1(I)*DFACT
368 GO TO 861
369862 CONTINUE
370831 CONTINUE
371832 CONTINUE
372 VACDST=VACDST*DUNITO/DUNIT
373 DUNITO=DUNIT
374 DO 871 JR=1,6
375 MD=MED(JR)
376 IF ((MD.GE.1).AND.(MD.LE.NMED)) THEN
377 ECUT(JR)=MAX(ECUT(JR),AE,AP+1.1*RM)
378 PCUT(JR)=MAX(PCUT(JR),AP)
379 IF((RHOR(JR).EQ.0.0))RHOR(JR)=RHO
380 END IF
381871 CONTINUE
382872 CONTINUE
383 IF (NMED.EQ.1) THEN
384 WRITE(KMPO,880)
385880 FORMAT(' EGS SUCCESSFULLY ''HATCHED'' FOR ONE MEDIUM.')
386 ELSE
387 WRITE(KMPO,890)NMED
388890 FORMAT(' EGS SUCCESSFULLY ''HATCHED'' FOR ',I5,' MEDIA.')
389 END IF
390 RETURN
391720 WRITE(KMPO,900)KMPI
392900 FORMAT(' END OF FILE ON UNIT ',I2,//, ' PROGRAM STOPPED IN HATCH '
393 *, 'BECAUSE THE'/ ' FOLLOWING NAMES WERE NOT RECOGNIZED:',/)
394 DO 911 IM=1,NMED
395 IF (LOK.NE.1) THEN
396 WRITE(KMPO,920)(MEDIA(I:I),I=1,LMDN)
397920 FORMAT(40X,'''',24A1,'''')
398 END IF
399911 CONTINUE
400912 CONTINUE
401 STOP
402 END
Note: See TracBrowser for help on using the repository browser.