source: branches/start/MagicSoft/Simulation/Corsika/Mmcs/howfar.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: 7.7 KB
Line 
1 SUBROUTINE HOWFAR
2C
3C*********************************************************************
4C DESIGN : D. HECK IK3 FZK KARLSRUHE
5C DATE : SEP 05, 1988
6C*********************************************************************
7C THE FOLLOWING IS A GENERAL SPECIFICATION OF HOWFAR:
8C GIVEN A PARTICLE AT (X,Y,Z) IN REGION IR AND GOING IN DIRECTION
9C (U,V,W), THIS ROUTINE ANSWERS THE QUESTION, CAN THE PARTICLE GO
10C A DISTANCE USTEP WITHOUT CROSSING A BOUNDARY OR OBSERVATION LEVEL?
11C IF YES, IT CALCULATES DNEAR AND RETURNS.
12C IF NO, IT SETS USTEP=DISTANCE TO BOUNDARY OR DETECTOR IN
13C IN THE CURRENT DIRECTION.
14C IT SETS IRNEW TO THE REGION NUMBER ON THE FAR SIDE
15C OF THE BOUNDARY (THIS CAN BE MESSY IN GENERAL!);
16C IT SETS NEWOBS TO THE DETECTOR NUMBER NEXT AFTER THE
17C DETECTOR JUST PASSING.
18C THE USER CAN TERMINATE A HISTORY BY SETTING IDISC>0. THE USER
19C CAN TRANSPORT THE LAST PARTICLE BY SETTING IDISC<0. HERE WE
20C TERMINATE ALL HISTORIES WHICH ENTER REGION 6 OR ARE GOING
21C BACKWARDS IN REGION 1 OR HAVE PASSED THE LAST OBSERVATION LEVEL.
22C*********************************************************************
23C ELECTRON OR PHOTON POSITIVE Z-DIRECTION (W>0) IS DOWNWARDS
24C |
25C | REGION 1 (VACUUM)
26C V
27C--------------------------- STARTING PLANE AT -BOUND(1) = -ZALTIT
28C
29C REGION 2 (AIR WITH EXPONENTIALLY
30C INCREASING DENSITY)
31C
32C--------------------------- BOUNDARY AT -BOUND(2)
33C
34C REGION 3 (AIR WITH EXPONENTIALLY
35C INCREASING DENSITY)
36C
37C--------------------------- BOUNDARY AT -BOUND(3)
38C
39C REGION 4 (AIR WITH EXPONENTIALLY
40C INCREASING DENSITY)
41C
42C--------------------------- BOUNDARY AT -BOUND(4)
43C
44C REGION 5 (AIR WITH EXPONENTIALLY
45C INCREASING DENSITY)
46C
47C-------------------------Z=0 BOUNDARY AT -BOUND(5) (SEA LEVEL)
48C////////////|/////////
49C////////////|///////// REGION 6 (VACUUM)
50C////////////V///////// (MAY CONTAIN DETECTOR)
51C ELECTRON OR PHOTON
52C------------------------ BOUNDARY AT -BOUND(6)
53C
54C*********************************************************************
55*KEEP,EPCONT.
56 COMMON/EPCONT/ EDEP,RATIO,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,IDISC,
57 * IROLD,IRNEW,RHOFAC, EOLD,ENEW,EKE,ELKE,BETA2,GLE,
58 * TSCAT,IAUSFL
59 DOUBLE PRECISION EDEP,RATIO
60 REAL TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,RHOFAC,EOLD,ENEW,
61 * EKE,ELKE,BETA2,GLE,TSCAT
62 INTEGER IDISC,IROLD,IRNEW,IAUSFL(29)
63*KEND.
64 COMMON/GEOM/ZALTIT,BOUND(6),NEWOBS,OBSLVL(10)
65*KEEP,OBSPAR.
66 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
67 * THETPR,PHIPR,NOBSLV
68 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
69 * THETAP,THETPR(2),PHIP,PHIPR(2)
70 INTEGER NOBSLV
71*KEEP,RUNPAR.
72 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
73 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
74 * MONIOU,MDEBUG,NUCNUC,
75 * CETAPE,
76 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
77 * N1STTR,MDBASE,
78 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
79 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
80 * ,GHEISH,GHESIG
81 COMMON /RUNPAC/ DSN,HOST,USER
82 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
83 REAL STEPFC
84 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
85 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
86 * N1STTR,MDBASE
87 INTEGER CETAPE
88 CHARACTER*79 DSN
89 CHARACTER*20 HOST,USER
90
91 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
92 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
93 * ,GHEISH,GHESIG
94*KEEP,STACKE.
95 COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
96 DOUBLE PRECISION E(60),TIME(60)
97 REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
98 INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
99*KEND.
100 COMMON/ACLOCK/NCLOCK,JCLOCK
101C_____IF (NCLOCK.GT.JCLOCK) THEN
102C______WRITE(MDEBUG,* )' HOWFAR:NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
103C______CALL AUSGB2
104C_____END IF
105 IF (IR(NP).GT.1 .AND. IR(NP).LT.6) THEN
106C *** WE ARE IN THE ATMOSPHERE - CHECK THE GEOMETRY
107 IRL=IR(NP)
108C *** GOING FORWARD - CONSIDER FIRST SINCE MOST FREQUENT
109 NOBS=IOBS(NP)
110 IF (W(NP).GT.0.0) THEN
111C *** TVAL IS DISTANCE TO NEXT BOUNDARY OR
112C *** OBSERVATION LEVEL IN THIS DIRECTION
113 TVAL=(-Z(NP)-MAX(BOUND(IRL),OBSLVL(NOBS)))/W(NP)
114 IF (TVAL.GT.USTEP) THEN
115C *** CAN TAKE CURRENTLY REQUESTED STEP
116 DNEAR(NP)=TVAL*W(NP)
117 ELSE
118C *** GO TO DETECTOR OR BOUNDARY, WHICH IS CLOSER
119 USTEP=MAX(TVAL,0.0001)
120 IF (BOUND(IRL).GE.OBSLVL(NOBS)) THEN
121C *** PARTICLE CROSSES BOUNDARY
122 IRNEW=IRL+1
123C *** PARTICLE LEAVES AIR
124 IF((IRNEW.GE.6))IDISC=-1
125 END IF
126 IF (BOUND(IRL).LE.OBSLVL(NOBS)) THEN
127C *** PARTICLE CROSSES DETECTOR
128 NEWOBS=NOBS+1
129C *** MAKE A VERY SMALL STEP TO AVOID HANGUP OF PROGRAM
130 IF((USTEP.LE.0.0))USTEP = 0.0001
131C *** TRANSPORT PARTICLE TO FINAL DETECTOR LEVEL AND DISCARD IT
132 IF((NEWOBS.GT.NOBSLV))IDISC=-1
133 END IF
134 END IF
135C *** END OF W(NP)>0 CASE
136C *** GOING UPWARD IN ATMOSPHERE
137 ELSE IF(W(NP).LT.0.0) THEN
138C *** NO DETECTOR ABOVE PARTICLE
139 IF (NOBS.LE.1) THEN
140C *** DISTANCE TO BOUNDARY ABOVE
141 TVAL=(-Z(NP)-BOUND(IRL-1))/W(NP)
142 IF (TVAL.GT.USTEP) THEN
143C *** CAN TAKE CURRENTLY REQUESTED STEP
144 DNEAR(NP)=MIN(Z(NP)+BOUND(IRL-1),-(Z(NP)+BOUND(IRL)))
145 ELSE
146C *** CROSS BOUNDARY ABOVE
147 USTEP=MAX(TVAL,0.0001)
148 IRNEW=IRL-1
149 END IF
150 ELSE
151C *** BOUNDARY AND DETECTOR ABOVE PARTICLE
152 TVAL=(-Z(NP)-MIN(BOUND(IRL-1),OBSLVL(NOBS-1)))/W(NP)
153 IF (TVAL.GT.USTEP) THEN
154C *** CAN TAKE CURRENTLY REQUESTED STEP
155C *** DNEAR IS CLOSEST DISTANCE TO DETECTOR OR
156C *** BOUNDARY ABOVE OR BELOW PARTICLE
157 DNEAR(NP)=MIN(Z(NP)+MIN(BOUND(IRL-1),OBSLVL(NOBS-1)), -Z(NP) +
158 * MAX(BOUND(IRL),OBSLVL(NOBS)))
159 ELSE
160C *** TAKE ONLY STEP UP TO BOUNDARY OR DETECTOR
161 USTEP=MAX(TVAL,0.0001)
162 IF (BOUND(IRL-1).LE.OBSLVL(NOBS-1)) THEN
163C *** PARTICLE CROSSES BOUNDARY ABOVE
164 IRNEW=IRL-1
165C *** PARTICLE LEAVES ATMOSPHERE
166 IF((IRNEW.LE.1))IDISC=1
167 END IF
168 IF ((BOUND(IRL-1).GE.OBSLVL(NOBS-1))) THEN
169C *** PARTICLE CROSSES DETECTOR ABOVE; IT IS NOT
170C *** PRINTED, BECAUSE IT MUST HIT DETECTOR DOWNWARDS
171 NEWOBS=NOBS-1
172 IOBS(NP)=NEWOBS
173 END IF
174 END IF
175 END IF
176C *** END W(NP)<0 CASE
177C *** PARTICLE IS MOVING HORIZONTALLY, CANNOT HIT BOUNDARY
178 ELSE IF(W(NP).EQ.0.0) THEN
179 RETURN
180 END IF
181C *** END OF ATMOSPHERE REGION CASE
182 ELSE IF(IR(NP).EQ.6) THEN
183C *** TERMINATE THIS HISTORY, IT IS PAST THE ATMOSPHERE
184 IDISC=1
185C *** WE ARE IN THE REGION WITH SOURCE ABOVE AIR
186 ELSE IF(IR(NP).EQ.1) THEN
187 IF (W(NP).GT.0.0) THEN
188C *** IT MUST BE A SOURCE PARTICLE ON BOUNDARY 1
189 USTEP=0.0001
190 IRNEW=2
191 ELSE
192C *** IT IS A REFLECTED PARTICLE, DISCARD IT
193 IDISC=1
194 END IF
195C *** END REGION 1 CASE
196 END IF
197 RETURN
198 END
Note: See TracBrowser for help on using the repository browser.