]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/peanut/prepre.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Minor corrections after big transformer changes
[u/mrichter/AliRoot.git] / GEANT321 / peanut / prepre.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1995/10/24 10:22:03 cernlib
6* Geant
7*
8*
9#include "geant321/pilot.h"
10*CMZ : 3.21/02 28/03/94 01.30.59 by S.Giani
11*-- Author :
12#if defined(CERNLIB_HPUX)
13$OPTIMIZE OFF
14#endif
15*$ CREATE PREPRE.FOR
16*COPY PREPRE
17*
18*=== prepre ===========================================================*
19*
20 SUBROUTINE PREPRE ( WEE, NNEUT, NPROT, NHOLE, ICYCL )
21
22#include "geant321/dblprc.inc"
23#include "geant321/dimpar.inc"
24#include "geant321/iounit.inc"
25*
26*----------------------------------------------------------------------*
27*----------------------------------------------------------------------*
28*
29#include "geant321/balanc.inc"
30#include "geant321/eva0.inc"
31#include "geant321/fheavy.inc"
32#include "geant321/finuc.inc"
33#include "geant321/nucdat.inc"
34#include "geant321/nucgeo.inc"
35#include "geant321/parevt.inc"
36#include "geant321/parnuc.inc"
37#include "geant321/part.inc"
38#include "geant321/resnuc.inc"
39*
40 REAL RNDM(2)
41 COMMON / FKCOSP / C1ST (3), C2ND (3), LEMISS
42 LOGICAL LEMISS
43 COMMON / FKCMCY / NPCYCL (20,2), IEVT, LOUT
44 COMMON / FKPLOC / IABCOU
45*
46 DIMENSION ERNCM (2), EPSMX (2), AMMAFT (2), AMNAFT (2),
47 & ZAFTR (2), EEXMNM (2), EEXDLT (2), EEXANW (2),
48 & EPSEEX (2), EPSDLT (2), EPSANY (2), EPSFIX (29),
49 & C (3), UMULEG (0:2), ACOLEG (0:2), NPART (2),
50 & NPRFIX (2)
51 LOGICAL LSPREJ, LNUCTS
52*
53 NPART (1) = NNEUT
54 NPART (2) = NPROT
55 IF ( LHLFIX ) THEN
56 NHLFIX = NHLEXP
57 NHOLE = NHOLE - NHLFIX
58 EHLFIX = 0.D+00
59 DO 50 IHOLE = 1, NHLFIX
60 EHLFIX = EHLFIX + HOLEXP (IHOLE)
61 50 CONTINUE
62 IF ( NHOLE .GT. 0 ) THEN
63 RHOIMP = ( RHOIMP * ( NHOLE + NHLFIX ) - RHOEXP ) / NHOLE
64 EKFIMP = ( EKFIMP * ( NHOLE + NHLFIX ) - EKFEXP ) / NHOLE
65 ELSE
66 RHOIMP = RHOAVE
67 EKFIMP = EKFAVE
68 END IF
69 ELSE
70 NHLFIX = 0
71 EHLFIX = 0.D+00
72 END IF
73 NHINI = NHOLE
74 IF ( PTRES .LT. ANGLGB ) THEN
75 UMO2 = ERES * ERES
76 UMO = SQRT (UMO2)
77 EKRES = 0.D+00
78 GAMCM = 1.D+00
79 ETAX = 0.D+00
80 ETAY = 0.D+00
81 ETAZ = 0.D+00
82 PXORI = 0.D+00
83 PYORI = 0.D+00
84 PZORI = 0.D+00
85 PCMORI = 0.D+00
86 CALL RACO ( CXAXCM, CYAXCM, CZAXCM )
87 ELSE
88 UMO2 = ( ERES - PTRES ) * ( ERES + PTRES )
89 UMO = SQRT (UMO2)
90 EKRES = ERES - UMO
91 GAMCM = ERES / UMO
92 ETACM = PTRES / UMO
93 ETAX = PXRES / UMO
94 ETAY = PYRES / UMO
95 ETAZ = PZRES / UMO
96 ETAPCM = ETAX * PXORI + ETAY * PYORI + ETAZ * PZORI
97 PHELP = EKORI + AM (KPORI) - ETAPCM / ( GAMCM + 1.D+00 )
98 PCMSX = PXORI - ETAX * PHELP
99 PCMSY = PYORI - ETAY * PHELP
100 PCMSZ = PZORI - ETAZ * PHELP
101 PCMORI = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
102 CXAXCM = PCMSX / PCMORI
103 CYAXCM = PCMSY / PCMORI
104 CZAXCM = PCMSZ / PCMORI
105 END IF
106 BBTAR = IBTAR
107 ZZTAR = ICHTAR
108 TVCMS = UMO - AMNRES
109 EHLFIX = MIN ( TVCMS, EHLFIX )
110 TVEFF = TVCMS - EHLFIX
111 IF ( TVEFF .LE. 0.D+00 ) GO TO 7000
112 TEMNU = SQRT ( TVCMS / ANOW * ALEVEL )
113 NEXMX = SQRT ( 2.D+00 * ANOW * TVCMS / ALEVEL )
114 NPTOT = NPART (1) + NPART (2)
115 IF ( NHOLE + NPTOT + NHLFIX .GE. NEXMX .OR.
116 & NPTOT .GT. NINT (0.5D+00 * ANOW) ) GO TO 7000
117 ANPTOT = NPTOT
118 ANPROT = NPART (2)
119 ANNEUT = NPART (1)
120 AVEBIN = ( ( BBTAR - ZZTAR - ACOLL + ZCOLL ) * BNENRG (2)
121 & + ( ZZTAR - ZCOLL ) * BNENRG (1) ) / ( BBTAR - ACOLL )
122 IAAFT = IBRES-1
123 IZAFT = ICRES
124 AAFTR = IAAFT
125 ZAFTR (1) = IZAFT
126 AMMAFT(1) = AAFTR * AMUAMU + 1.D-03 * FKENER (AAFTR,ZAFTR(1))
127 AMNAFT(1) = AMMAFT (1) - ZAFTR (1) * AMELEC + ELBNDE (IZAFT)
128 CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (1), EEXMNM (1), EEXDUM )
129 BNDGAV = BNENRG (2)
130 EPSMX (1) = UMO - AMNRES - BNDGAV - EHLFIX
131 EPSDLT (1) = UMO - AMNRES - BNDGAV - EEXDLT (1)
132 EPSEEX (1) = UMO - AMNRES - BNDGAV - EEXMNM (1)
133 IZAFT = ICRES - 1
134 ZAFTR (2) = IZAFT
135 AMMAFT(2) = AAFTR * AMUAMU + 1.D-03 * FKENER (AAFTR,ZAFTR(2))
136 AMNAFT(2) = AMMAFT (2) - ZAFTR (2) * AMELEC + ELBNDE (IZAFT)
137 CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (2), EEXMNM (2), EEXDUM )
138 BNDGAV = BNENRG (1)
139 EPSMX (2) = UMO - AMNRES - BNDGAV - EHLFIX
140 IF ( EPSMX (1) + EPSMX (2) .LT. 2.D+00 * TEMNU ) GO TO 7000
141 EPSDLT (2) = UMO - AMNRES - BNDGAV - EEXDLT (2)
142 EPSEEX (2) = UMO - AMNRES - BNDGAV - EEXMNM (2)
143 IF ( NP .LE. NP0 ) THEN
144 IF ( KPORI .EQ. 8 ) THEN
145 EEXANW (1) = EEXDLT (1)
146 EPSANY (1) = MIN ( EPSDLT (1), EPSMX (1) )
147 EEXANW (2) = 0.D+00
148 EPSANY (2) = EPSMX (2)
149 ELSE IF ( KPORI .EQ. 1 ) THEN
150 EEXANW (1) = 0.D+00
151 EPSANY (1) = EPSMX (1)
152 EEXANW (2) = EEXDLT (2)
153 EPSANY (2) = MIN ( EPSDLT (2), EPSMX (2) )
154 ELSE
155 EEXANW (1) = 0.D+00
156 EPSANY (1) = EPSMX (1)
157 EEXANW (2) = 0.D+00
158 EPSANY (2) = EPSMX (2)
159 END IF
160 ELSE
161 EEXANW (1) = 0.D+00
162 EPSANY (1) = EPSMX (1)
163 EEXANW (2) = 0.D+00
164 EPSANY (2) = EPSMX (2)
165 END IF
166 SIGIN0 = PI * ( R0SIGK * RMASS (IBRES-1) )**2
167 NEMISS = 0
168 MNUCTS = 0
169 1000 CONTINUE
170 NPRFIX (1) = 0
171 NPRFIX (2) = 0
172 JEMFIX = 0
173 EPRFIX = 0.D+00
174 IF ( .NOT. LEMISS ) ICYCL = ICYCL + 1
175 IF ( NPTOT .LE. 0 ) GO TO 4600
176 IF ( NNUCTS .GT. MNUCTS ) THEN
177 MNUCTS = MNUCTS + 1
178 JNUCTS = INUCTS (MNUCTS)
179 JEMIS0 = 2 - ABS (KPNUCL(JNUCTS)) / 8
180 JEMIS1 = JEMIS0
181 JDEMIS = 1
182 EKFSAV = EKFIMP
183 RHOSAV = RHOIMP
184 TMPRHO = 0.001D+00 * RHOAVE
185 RHOIMP = MAX ( RHNUCL (JNUCTS), TMPRHO )
186 TMPEKF = 0.001D+00 * EKFAVE
187 EKFIMP = MAX ( EKFNUC (JNUCTS), TMPEKF )
188 LNUCTS = .TRUE.
189 ELSE
190 LNUCTS = .FALSE.
191 CALL GRNDM(RNDM,1)
192 IF ( RNDM (1) .LT. ANNEUT / ANPTOT ) THEN
193 JEMIS0 = 1
194 JEMIS1 = 2
195 JDEMIS = 1
196 ELSE
197 JEMIS0 = 2
198 JEMIS1 = 1
199 JDEMIS = -1
200 END IF
201 END IF
202 NEXTOT = NPTOT + NHOLE
203 RHORED = ACOLL / BBTAR
204 IF ( LGDHPR ) THEN
205 EKFCNN = ( ZNOW * EKFCEN (1) + ( ANOW - ZNOW ) * EKFCEN (2)
206 & ) / ANOW
207 RHOWEI = 0.5D+00
208 IF ( LNUCTS ) THEN
209 AEFREF = EKFIMP
210 AEFRAV = EKFIMP
211 RHONUC = RHOIMP
212 ELSE IF ( NPTOT .LE. 2 .AND. ICYCL .EQ. 1 ) THEN
213 AHLTOT = MAX ( NHOLE, 1 )
214 WEIGH1 = NHINI / AHLTOT
215 WEIGH2 = 1.D+00 - WEIGH1
216 AEFREF = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
217 RHONUC = WEIGH1 * RHOIMP + WEIGH2 * RHOAVE
218 AEFRAV = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
219 RHONUC = RHORED * ( ( AHLTOT + 1.D+00 - RHOWEI )
220 & * RHONUC + RHOWEI * RHOAVE ) / ( AHLTOT + 1.D+00 )
221 AEFRAV = ( AHLTOT * AEFRAV + EKFAVE )
222 & / ( AHLTOT + 1.D+00 )
223 ELSE
224 AHLTOT = MAX ( NHOLE, 1 )
225 WEIGH1 = NHINI / AHLTOT
226 WEIGH2 = 1.D+00 - WEIGH1
227 AEFREF = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
228 RHONUC = WEIGH1 * RHOIMP + WEIGH2 * RHOAVE
229 AEFRAV = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
230 RHONUC = RHORED * ( ( AHLTOT + 1.D+00 - RHOWEI )
231 & * RHONUC + RHOWEI * RHOAVE ) / ( AHLTOT + 1.D+00 )
232 AEFRAV = ( AHLTOT * AEFRAV + EKFAVE )
233 & / ( AHLTOT + 1.D+00 )
234 END IF
235 CLAMDI = 1.0D+00
236 ELSE
237 EKFCNN = AEFRMX
238 AEFREF = 0.5D+00 * AEFRMX
239 AEFRAV = AEFREF
240 RHONUC = 0.5D+00**1.5D+00 * RHONU0
241 CLAMDI = 0.5D+00
242 END IF
243 DO 4500 JEMISS = JEMIS0, JEMIS1, JDEMIS
244 IF ( NPART (JEMISS) .LE. 0 ) GO TO 4500
245 IEMISS = 3 - JEMISS
246 RATEC0 = 2.D+00 * AMNUCL (IEMISS) / PLABRC**3 / PISQ
247 BNDGEN = AMNAFT (JEMISS) + AMNUCL (IEMISS) - AMNRES
248 IF ( .NOT. LNUCTS .AND. EPRFIX .GT. ANGLGB ) THEN
249 TVEFF = TVEFF - EPRFIX
250 NPTOT = NPTOT - NPRFIX (1) - NPRFIX (2)
251 NPART (1) = NPART (1) - NPRFIX (1)
252 NPART (2) = NPART (2) - NPRFIX (2)
253 JEMFIX = JEMISS
254 ERNCM0 = ERNCM (JEMISS)
255 EPSMX0 = EPSMX (JEMISS)
256 BNDGAV = BNENRG (IEMISS)
257 EPSMX (JEMISS) = UMO - AMNRES - BNDGAV - EHLFIX - EPRFIX
258 ELSE
259 JEMFIX = 0
260 END IF
261 BNDGAV = BNENRG (IEMISS)
262 EPSMIN = BNDGEN - BNDGAV
263 BNDHLP = BNENRG (IEMISS)
264 IF ( EPSMX (JEMISS) - EPSMIN .LT. TEMNU ) GO TO 4500
265 IF ( EHLFIX + EPRFIX .LE. 0.5D+00 * EEXANW (JEMISS) ) THEN
266 UUMIN = TVEFF - 0.5D+00 * ( EPSANY (JEMISS)
267 & + EPSMX (JEMISS) ) - BNDGAV
268 ELSE
269 UUMIN = TVEFF - EPSMX (JEMISS) - BNDGAV
270 END IF
271 UUMAX = TVEFF - BNDGEN
272 BNDUSE = BNENRG (IEMISS)
273 IF ( LNUCTS ) THEN
274 EKNNUC = ENNUC (JNUCTS) - AMNUCL (IEMISS) + BNDGEN
275 & - RSTNUC (JNUCTS)
276 EPSHLP = BNDHLP - BNDGEN + EPSMIN
277 EKNNUC = EKNNUC + EPSHLP
278 IF ( EKNNUC .LE. EPSMIN ) THEN
279 GO TO 4500
280 ELSE IF ( EKNNUC .LE. 0.D+00 ) THEN
281 DELTAE = 0.5D+00 * ( - EPSMIN - EKNNUC )
282 EKNNUC = EKNNUC + DELTAE
283 ELSE
284 DELTAE = 0.D+00
285 END IF
286 ENNUC (JNUCTS) = EKNNUC + AMNUCL (IEMISS)
287 PNUCCO = SQRT ( EKNNUC * ( EKNNUC + 2.D+00
288 & * AMNUCL (IEMISS) ) )
289 PHELP = PNUCCO / PNUCL (JNUCTS)
290 PCMSX = PXNUCL (JNUCTS) * PHELP
291 PCMSY = PYNUCL (JNUCTS) * PHELP
292 PCMSZ = PZNUCL (JNUCTS) * PHELP
293 ETAPCM = ETAX * PCMSX + ETAY * PCMSY + ETAZ * PCMSZ
294 PHELP = ENNUC (JNUCTS) - ETAPCM / ( GAMCM + 1.D+00 )
295 PCMSX = PCMSX - ETAX * PHELP
296 PCMSY = PCMSY - ETAY * PHELP
297 PCMSZ = PCMSZ - ETAZ * PHELP
298 PCMST = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
299 EPS = GAMCM * ENNUC (JNUCTS) - ETAPCM - AMNUCL (IEMISS)
300 IF ( EPS .LE. EPSMIN ) THEN
301 GO TO 4500
302 ELSE IF ( EPS .LT. ANGLGB * ENNUC (JNUCTS) ) THEN
303 GO TO 4500
304 END IF
305 EPS = EPS - DELTAE
306 C (1) = PCMSX / PCMST
307 C (2) = PCMSY / PCMST
308 C (3) = PCMSZ / PCMST
309 EPS = MIN ( EPS, HLFHLF * ( EPSMX (JEMISS)
310 & + EPSANY (JEMISS) ) )
311 IF ( JEMISS .EQ. 2 ) THEN
312 FLKCOU = DOST ( 1, ZAFTR (JEMISS) )
313 CCOUL = DOST ( 3, ZAFTR (JEMISS) )
314 ETHRES = FLKCOU * COULBH * ( ZNOW - 1.D+00 )
315 & / RMASS (IBRES-1)
316 END IF
317 FLAMMX = 1.D+00
318 NTENT = 1
319 ISAMPL = 0
320 LSPREJ = .FALSE.
321 GO TO 3000
322 END IF
323 IF ( UUMIN .GE. UUMAX ) GO TO 4500
324 EPSINV = 0.5D+00 * ( UMO2 - ( AMNUCL (IEMISS)
325 & + AMNAFT (JEMISS) )**2 ) / AMNAFT (JEMISS)
326 EPSWLL = EPSINV + AEFRAV + BNDUSE + AMNUCL (IEMISS)
327 BETWLL = SQRT ( 1.D+00 - AMNUSQ (IEMISS) / EPSWLL**2 )
328 EKEWLL = EPSWLL - AMNUCL (IEMISS)
329 IF ( JEMISS .EQ. 1 ) THEN
330 AALPHA = 0.76D+00 + 2.2D+00 / RMASS (IBRES-1)
331 BBETAA = ( 2.22D+00 / RMASS (IBRES-1)**2 - 0.05D+00 )
332 & / AALPHA * 1.D-03
333 TMP07 = 0.07D+00
334 EPSEFF = MIN ( EPSINV, TMP07 )
335 SIGINV = 0.1D+00 * XINNEU ( EPSEFF, ZAFTR (JEMISS),
336 & BBETAA )
337 SFTYMX = 1.1D+00
338 GSNGLV = 1.5D+00 * ( ANOW - ZNOW ) / EKFCEN (IEMISS)
339 IF ( EKEWLL .LE. 0.04D+00 ) THEN
340 SIGMAP = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
341 & -1.86D+00 + 0.09415D+03 * EKEWLL
342 & + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
343 & * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
344 & + 0.13D+03 * EKEWLL )**2 )
345 IF ( EKEWLL .LT. 0.02D+00 ) THEN
346 SIGMAN = 0.3333333333333333D+00 * SIGMAP
347 ELSE
348 SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
349 & + 42.9D+00
350 END IF
351 ELSE
352 SIGMAP = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
353 & + 82.2D+00
354 SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
355 & + 42.9D+00
356 END IF
357 SIGMAP = 0.1D+00 * SIGMAP
358 SIGMAN = 0.1D+00 * SIGMAN
359 SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
360 & + ( ZNOW - ANPROT ) * SIGMAP )
361 & / ( ANOW - ANPTOT )
362 ELSE
363 EPSCOU = UMO - AMNAFT (JEMISS) - AMNUCL (IEMISS)
364 FLKCOU = DOST ( 1, ZAFTR (JEMISS) )
365 CCOUL = DOST ( 3, ZAFTR (JEMISS) )
366 ETHRES = FLKCOU * COULBH * ( ZNOW - 1.D+00 )
367 & / RMASS (IBRES-1)
368 TMP07 = 0.07D+00
369 EPSEFF = MIN ( EPSINV, TMP07 )
370 SIGINV = 1.D-01 * XINPRO ( EPSEFF, ZAFTR (JEMISS),
371 & ETHRES )
372 SFTYMX = 1.2D+00
373 GSNGLV = 1.5D+00 * ZNOW / EKFCEN (IEMISS)
374 IF ( EKEWLL .LE. 0.04D+00 ) THEN
375 SIGMAN = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
376 & -1.86D+00 + 0.09415D+03 * EKEWLL
377 & + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
378 & * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
379 & + 0.13D+03 * EKEWLL )**2 )
380 IF ( EKEWLL .LT. 0.02D+00 ) THEN
381 SIGMAP = 0.3333333333333333D+00 * SIGMAN
382 ELSE
383 SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
384 & + 42.9D+00
385 END IF
386 ELSE
387 SIGMAN = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
388 & + 82.2D+00
389 SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
390 & + 42.9D+00
391 END IF
392 SIGMAP = 0.1D+00 * SIGMAP
393 SIGMAN = 0.1D+00 * SIGMAN
394 SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
395 & + ( ZNOW - ANPROT ) * SIGMAP )
396 & / ( ANOW - ANPTOT )
397 END IF
398 ZITA = AEFRAV / EKEWLL
399 IF ( ZITA .LE. 0.5D+00 ) THEN
400 PZITA = 1.D+00 - 1.4D+00 * ZITA
401 ELSE
402 PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
403 & * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
404 END IF
405 RATEC = SIGINV * RATEC0 * EPSINV / GSNGLV
406 IF ( RATEC .LE. 0.D+00 ) GO TO 4500
407 ALAMDC = BETWLL / RATEC
408 ALAMDI = 1.D+00 / ( SGNUNU * RHONUC * PZITA ) + MAX (TWOTWO
409 & * PI * PLABRC * BETWLL * EPSWLL / AMNUSQ (IEMISS),
410 & RZNUCL )
411 AMUTOT = 1.D+00 / ALAMDI + 1.D+00 / ALAMDC
412 RATEI = BETWLL / ALAMDI
413 FLAMMX = SFTYMX * RATEC / ( RATEC + CLAMDI * RATEI )
414 IF ( FLAMMX .LT. 0.D+00 ) THEN
415 GO TO 4500
416 END IF
417 LSPREJ = .NOT. ( NHOLE .GT. 2 .OR. NPTOT .GT. 4 )
418 IF ( LSPREJ ) THEN
419 IF ( NPTOT .EQ. 1 ) THEN
420 IF ( NHOLE .EQ. 0 ) THEN
421 EPS = EPSMX (JEMISS)
422 FLAMMX = 1.D+00
423 NTENT = 1
424 ISAMPL = 0
425 LSPREJ = .FALSE.
426 GO TO 3000
427 ELSE IF ( NHOLE .EQ. 1 ) THEN
428 ISAMPL = 1
429 ELSE IF ( NHOLE .EQ. 2 ) THEN
430 ISAMPL = 2
431 ELSE
432 ISAMPL = 7
433 END IF
434 ELSE IF ( NPTOT .EQ. 2 ) THEN
435 IF ( NHOLE .EQ. 0 ) THEN
436 ISAMPL = 7
437 ELSE IF ( NHOLE .EQ. 1 ) THEN
438 ISAMPL = 3
439 ELSE IF ( NHOLE .EQ. 2 ) THEN
440 ISAMPL = 4
441 ELSE
442 ISAMPL = 7
443 END IF
444 ELSE IF ( NPTOT .EQ. 3 ) THEN
445 IF ( NHOLE .EQ. 0 ) THEN
446 ISAMPL = 7
447 ELSE IF ( NHOLE .EQ. 1 ) THEN
448 ISAMPL = 5
449 ELSE IF ( NHOLE .EQ. 2 ) THEN
450 ISAMPL = 6
451 ELSE
452 ISAMPL = 7
453 END IF
454 ELSE
455 ISAMPL = 7
456 END IF
457 ELSE
458 ISAMPL = 7
459 END IF
460 GO TO ( 2100, 2200, 2300, 2400, 2500, 2600, 2700 ), ISAMPL
461 2100 CONTINUE
462 LSPREJ = .FALSE.
463 AEFRPA = AEFREF
464 UUMAX = MIN ( AEFRPA, UUMAX )
465 AITOT = ( AEFRPA**3 - ( AEFRPA + UUMIN - UUMAX )**3 )
466 & / AEFRPA**3
467 IF ( UUMAX .GT. TVEFF ) THEN
468 AIPRO = MIN ( AITOT, ONEONE ) * FLAMMX
469 AITOT = AITOT * FLAMMX
470 ELSE
471 AITOT = AITOT * FLAMMX
472 AIPRO = AITOT
473 END IF
474 GO TO 2800
475 2200 CONTINUE
476 LSPREJ = .FALSE.
477 AEFRPA = AEFREF
478 UUMAX = MIN ( TWOTWO * AEFRPA, UUMAX )
479 IF ( TVEFF .LE. AEFRPA ) THEN
480 EDENOM = TVEFF**2 * 0.5D+00
481 ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
482 EDENOM = ( TVEFF**2 - 2.D+00 * ( TVEFF - AEFRPA )**2 )
483 & * 0.5D+00
484 ELSE
485 EDENOM = AEFRPA**2
486 END IF
487 UUMN2 = UUMIN**2
488 IF ( UUMAX .GT. AEFRPA ) THEN
489 UUDIV = AEFRPA
490 UUDV2 = UUDIV**2
491 AIPR1 = 0.5D+00 * UUDV2
492 AIPR2 = 0.5D+00 * ( UUMAX - UUDIV ) * ( 3.D+00 * UUDIV
493 & - UUMAX )
494 IF ( UUMIN .LE. UUDIV ) THEN
495 AITO1 = 0.5D+00 * ( UUDV2 - UUMN2 )
496 AITO2 = AIPR2
497 AIHLP = 0.D+00
498 ELSE
499 AITO1 = 0.D+00
500 AIHLP = 0.5D+00 * ( UUMIN - UUDIV ) * ( 3.D+00
501 & * UUDIV - UUMIN )
502 AITO2 = AIPR2 - AIHLP
503 END IF
504 ELSE
505 UUDIV = UUMAX
506 UUDV2 = UUDIV**2
507 AIPR1 = 0.5D+00 * UUDV2
508 AITO1 = 0.5D+00 * ( UUDV2 - UUMN2 )
509 AITO2 = 0.D+00
510 AIPR2 = 0.D+00
511 END IF
512 AITOT = AITO1 + AITO2
513 AITO1 = AITO1 / AITOT
514 AITO2 = AITO2 / AITOT
515 AITOT = AITOT * FLAMMX / EDENOM
516 IF ( UUMAX .GT. TVEFF ) THEN
517 AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM, ONEONE )
518 & * FLAMMX
519 ELSE
520 AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
521 END IF
522 GO TO 2800
523 2300 CONTINUE
524 LSPREJ = .FALSE.
525 AEFRPA = AEFREF
526 IF ( TVEFF .LE. AEFRPA ) THEN
527 EDENOM = 0.25D+00 * TVEFF * TVEFF
528 ELSE
529 EDENOM = 0.25D+00 * AEFRPA * ( 2.D+00 * TVEFF
530 & - AEFRPA )
531 END IF
532 FCHLP = 0.25D+00 * NPART (JEMISS)
533 IF ( UUMAX .LE. AEFRPA ) THEN
534 UUDIV = UUMAX
535 AIPR1 = FCHLP * UUDIV**2
536 AITO1 = AIPR1 - FCHLP * UUMIN**2
537 AITO2 = 0.D+00
538 AIPR2 = 0.D+00
539 ELSE
540 UUDIV = AEFRPA
541 AIPR1 = FCHLP * UUDIV**2
542 AIPR2 = FCHLP * ( UUMAX - UUDIV ) * UUDIV
543 IF ( UUMIN .GT. UUDIV ) THEN
544 AITO1 = 0.D+00
545 AIHLP = FCHLP * ( UUMIN - UUDIV ) * UUDIV
546 AITO2 = AIPR2 - AIHLP
547 ELSE
548 AIHLP = 0.D+00
549 AITO1 = AIPR1 - FCHLP * UUMIN**2
550 AITO2 = AIPR2
551 END IF
552 END IF
553 AITOT = AITO1 + AITO2
554 AITO1 = AITO1 / AITOT
555 AITO2 = AITO2 / AITOT
556 IF ( UUMAX .GT. TVEFF ) THEN
557 DDNPAR = NPART(JEMISS)
558 AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM,
559 & DDNPAR ) * FLAMMX
560 ELSE
561 AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
562 END IF
563 GO TO 2800
564 2400 CONTINUE
565 LSPREJ = .TRUE.
566 AEFRPA = AEFREF
567 IF ( TVEFF .LE. AEFRPA ) THEN
568 LSPREJ = .FALSE.
569 FEMAX = 1.D+00
570 ISAMPL = 7
571 ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
572 FEMAX = TVEFF**3 / ( TVEFF**3 - 2.D+00 * ( TVEFF
573 & - AEFRPA )**3 )
574 ELSE
575 FEMAX = 0.1666666666666667D+00 * TVEFF**3
576 & / ( AEFRPA**2 * ( TVEFF - AEFRPA ) )
577 END IF
578 UUMNR = ( UUMIN / TVEFF )**(NEXTOT-1)
579 UUMXR = ( UUMAX / TVEFF )**(NEXTOT-1)
580 AITOT = NPART (JEMISS) * FLAMMX * ( UUMXR - UUMNR )
581 & * FEMAX
582 AIPRO = NPART (JEMISS) * FLAMMX * UUMXR * FEMAX
583 GO TO 2800
584 2500 CONTINUE
585 LSPREJ = .FALSE.
586 AEFRPA = AEFREF
587 IF ( TVEFF .LE. AEFRPA ) THEN
588 EDENOM = TVEFF**3 / 36.D+00
589 ELSE
590 EDENOM = ( TVEFF**3 - ( TVEFF - AEFRPA )**3 )
591 & / 26.D+00
592 END IF
593 IF ( UUMAX .GT. AEFRPA ) THEN
594 UUDIV = AEFRPA
595 UUDV3 = UUDIV**3
596 UUMN3 = UUMIN**3
597 FCHLP = NPART(JEMISS) / 36.D+00
598 AIPR1 = FCHLP * UUDV3
599 AIPR2 = 3.D+00 * FCHLP * UUMAX * UUDIV * ( UUMAX
600 & - UUDIV )
601 IF ( UUMIN .LE. UUDIV ) THEN
602 AITO1 = AIPR1 - FCHLP * UUMN3
603 AITO2 = AIPR2
604 AIHLP = 0.D+00
605 ELSE
606 AITO1 = 0.D+00
607 AIHLP = 3.D+00 * FCHLP * UUMIN * UUDIV * ( UUMIN
608 & - UUDIV )
609 AITO2 = AIPR2 - AIHLP
610 END IF
611 ELSE
612 UUDIV = UUMAX
613 UUDV3 = UUDIV**3
614 UUMN3 = UUMIN**3
615 FCHLP = NPART(JEMISS) / 36.D+00
616 AIPR1 = FCHLP * UUDV3
617 AITO1 = AIPR1 - FCHLP * UUMN3
618 AIPR2 = 0.D+00
619 AITO2 = 0.D+00
620 END IF
621 AITOT = AITO1 + AITO2
622 AITO1 = AITO1 / AITOT
623 AITO2 = AITO2 / AITOT
624 AITOT = AITOT * FLAMMX / EDENOM
625 IF ( UUMAX .GT. TVEFF ) THEN
626 DDNPAR = NPART(JEMISS)
627 AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM,
628 & DDNPAR ) * FLAMMX
629 ELSE
630 AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
631 END IF
632 GO TO 2800
633 2600 CONTINUE
634 AEFRPA = AEFREF
635 IF ( TVEFF .LE. AEFRPA ) THEN
636 EDENOM = TVEFF**4 / 288.D+00
637 ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
638 EDENOM = ( TVEFF**4 - 2.D+00 * ( TVEFF - AEFRPA )**4 )
639 & / 288.D+00
640 ELSE
641 EDENOM = AEFRPA**2 * ( 0.5D+00 * TVEFF**2 - TVEFF
642 & * AEFRPA + 7.D+00 / 12.D+00 * AEFRPA**2 )
643 & / 12.D+00
644 END IF
645 FCHLP = NPART (JEMISS) / 288.D+00
646 IF ( UUMAX .GT. 2.D+00 * AEFRPA ) THEN
647 LSPREJ = .TRUE.
648 UUDIV = AEFRPA
649 UUDI2 = 2.D+00 * AEFRPA
650 UUDV4 = UUDIV**4
651 UUMN4 = UUMIN**4
652 UUD24 = UUDI2**4
653 AIPR1 = FCHLP * UUDV4
654 AIPR2 = FCHLP * ( UUD24 - UUDV4 )
655 AIPR3 = 12.D+00 * FCHLP * AEFRPA**2 * UUMAX
656 & * ( UUMAX - UUDI2 )
657 IF ( UUMIN .LE. UUDIV ) THEN
658 AITO1 = AIPR1 - FCHLP * UUMN4
659 AITO2 = AIPR2
660 AITO3 = AIPR3
661 AIHLP = 0.D+00
662 AIHL2 = 0.D+00
663 ELSE IF ( UUMIN .LE. UUDI2 ) THEN
664 AITO1 = 0.D+00
665 AIHLP = FCHLP * ( UUMN4 - UUDV4 )
666 AITO2 = AIPR2 - AIHLP
667 AITO3 = AIPR3
668 AIHL2 = 0.D+00
669 ELSE
670 AITO1 = 0.D+00
671 AITO2 = 0.D+00
672 AIHL2 = 12.D+00 * FCHLP * AEFRPA**2 * UUMIN
673 & * ( UUMIN - UUDI2 )
674 AITO3 = AIPR3 - AIHL2
675 END IF
676 ELSE IF ( UUMAX .GT. AEFRPA ) THEN
677 LSPREJ = .TRUE.
678 UUDIV = AEFRPA
679 UUDI2 = UUMAX
680 UUDV4 = UUDIV**4
681 UUMN4 = UUMIN**4
682 UUD24 = UUDI2**4
683 AIPR1 = FCHLP * UUDV4
684 AIPR2 = FCHLP * ( UUD24 - UUDV4 )
685 IF ( UUMIN .LE. UUDIV ) THEN
686 AITO1 = AIPR1 - FCHLP * UUMN4
687 AITO2 = AIPR2
688 AIHLP = 0.D+00
689 ELSE
690 AITO1 = 0.D+00
691 AIHLP = FCHLP * ( UUMN4 - UUDV4 )
692 AITO2 = AIPR2 - AIHLP
693 END IF
694 AIPR3 = 0.D+00
695 AITO3 = 0.D+00
696 ELSE
697 LSPREJ = .FALSE.
698 UUDI2 = UUMAX
699 UUDIV = UUMAX
700 UUDV4 = UUDIV**4
701 UUMN4 = UUMIN**4
702 AIPR1 = FCHLP * UUDV4
703 AITO1 = AIPR1 - FCHLP * UUMN4
704 AITO2 = 0.D+00
705 AIPR2 = 0.D+00
706 AITO3 = 0.D+00
707 AIPR3 = 0.D+00
708 END IF
709 AITOT = AITO1 + AITO2 + AITO3
710 AITO1 = AITO1 / AITOT
711 AITO2 = AITO2 / AITOT
712 AITO3 = AITO3 / AITOT
713 AITOT = AITOT * FLAMMX / EDENOM
714 AIPRO = ( AIPR1 + AIPR2 + AIPR3 ) * FLAMMX / EDENOM
715 GO TO 2800
716 2700 CONTINUE
717 LSPREJ = .FALSE.
718 TVLEFF = LOG(TVEFF)
719 IF(UUMAX.LE.0.) THEN
720 UULMXR=-100
721 ELSE
722 UULMXR = (LOG(UUMAX)-TVLEFF)*(NEXTOT-1)
723 ENDIF
724 IF(UULMXR.LT.-88) THEN
725 UUMNR = 0.
726 UUMXR = 0.
727 ELSE
728 UUMXR = EXP(UULMXR)
729 IF(UUMIN.LE.0.) THEN
730 UULMNR=-100
731 ELSE
732 UULMNR = (LOG(UUMIN)-TVLEFF)*(NEXTOT-1)
733 ENDIF
734 IF(UULMNR.LT.-88) THEN
735 UUMNR = 0.
736 ELSE
737 UUMNR = EXP(UULMNR)
738 ENDIF
739 ENDIF
740 AITOT = NPART (JEMISS) * FLAMMX * ( UUMXR - UUMNR )
741 AIPRO = NPART (JEMISS) * FLAMMX * MIN ( UUMXR, ONEONE )
742 2800 CONTINUE
743 NTENT = INT (AIPRO)
744 CALL GRNDM(RNDM,1)
745 RNTENT = RNDM (1)
746 IF ( RNTENT .LT. AIPRO - NTENT ) NTENT = NTENT + 1
747 3000 CONTINUE
748 ITACC = 0
749 DO 4100 IT = 1, NTENT
750 GO TO ( 3100, 3200, 3300, 3400, 3500, 3600, 3700 ),
751 & ISAMPL
752 GO TO 4000
753 3100 CONTINUE
754 CALL GRNDM(RNDM,1)
755 UURND = AEFRPA * ( 1.D+00 - AITOT / FLAMMX
756 & * RNDM (1) )**0.3333333333333333D+00
757 UURND = AEFRPA - UURND + UUMIN
758 GO TO 3800
759 3200 CONTINUE
760 CALL GRNDM(RNDM,1)
761 RNDEPS = RNDM (1)
762 IF ( RNDEPS .LT. AITO1 ) THEN
763 RNDEPS = RNDEPS / AITO1
764 UURND = SQRT ( RNDEPS * ( UUDV2 - UUMN2 ) + UUMN2)
765 ELSE
766 RNDEPS = RNDEPS - AITO1
767 RNDEPS = AIHLP + RNDEPS * AITOT * EDENOM / FLAMMX
768 UURND = 2.D+00 * UUDIV - SQRT ( UUDIV**2
769 & - 2.D+00 * RNDEPS )
770 END IF
771 GO TO 3800
772 3300 CONTINUE
773 CALL GRNDM(RNDM,1)
774 RNDEPS = RNDM (1)
775 IF ( RNDEPS .LT. AITO1 ) THEN
776 RNDEPS = RNDEPS * AITOT
777 UURND = RNDEPS / FCHLP + UUMIN**2
778 UURND = SQRT (UURND)
779 ELSE
780 RNDEPS = ( RNDEPS - AITO1 ) * AITOT + AIHLP
781 UURND = RNDEPS / FCHLP / UUDIV + UUDIV
782 END IF
783 GO TO 3800
784 3400 CONTINUE
785 CALL GRNDM(RNDM,1)
786 UURND = TVEFF * ( RNDM (1) * AITOT / FLAMMX
787 & / FEMAX / NPART(JEMISS) + UUMNR )**
788 & ( 1.D+00 / (NEXTOT-1) )
789 IF ( UURND .LE. AEFRPA ) THEN
790 LSPREJ = .FALSE.
791 FREJE = 1.D+00
792 ELSE IF ( UURND .LE. 2.D+00 * AEFRPA ) THEN
793 FREJE = ( UURND**2 - 2.D+00 * ( UURND - AEFRPA )
794 & **2 ) / UURND**2
795 ELSE
796 FREJE = 2.D+00 * AEFRPA**2 / UURND**2
797 END IF
798 GO TO 3800
799 3500 CONTINUE
800 CALL GRNDM(RNDM,1)
801 RNDEPS = RNDM (1)
802 IF ( RNDEPS .LT. AITO1 ) THEN
803 RNDEPS = RNDEPS / AITO1
804 UURND = RNDEPS * ( UUDV3 - UUMN3 ) + UUMN3
805 UURND = UURND**0.333333333333333D+00
806 LSPREJ = .FALSE.
807 ELSE IF ( RNDEPS .LT. AITO1 + AITO2 ) THEN
808 RNDEPS = RNDEPS - AITO1
809 RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHLP
810 UURND = 0.5D+00 * ( UUDIV + SQRT ( UUDIV**2
811 & + 1.333333333333333D+00 * RNDEPS / FCHLP
812 & / UUDIV ) )
813 LSPREJ = .FALSE.
814 END IF
815 GO TO 3800
816 3600 CONTINUE
817 CALL GRNDM(RNDM,1)
818 RNDEPS = RNDM (1)
819 IF ( RNDEPS .LT. AITO1 ) THEN
820 RNDEPS = RNDEPS / AITO1
821 UURND = RNDEPS * ( UUDV4 - UUMN4 ) + UUMN4
822 UURND = UURND**0.25D+00
823 LSPREJ = .FALSE.
824 ELSE IF ( RNDEPS .LT. AITO1 + AITO2 ) THEN
825 RNDEPS = RNDEPS - AITO1
826 RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHLP
827 UURND = RNDEPS / FCHLP + UUDV4
828 UURND = UURND**0.25D+00
829 LSPREJ = .TRUE.
830 FREJE = 1.D+00 - 2.D+00 * ( 1.D+00 - UUDIV / UURND
831 & )**3
832 ELSE
833 RNDEPS = RNDEPS - AITO1 - AITO2
834 RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHL2
835 UURND = AEFRPA + SQRT ( AEFRPA**2 + RNDEPS /
836 & ( 12.D+00 * FCHLP * AEFRPA**2 ) )
837 LSPREJ = .FALSE.
838 END IF
839 GO TO 3800
840 3700 CONTINUE
841 CALL GRNDM(RNDM,1)
842 UURND = TVEFF * ( RNDM (1) * AITOT / FLAMMX
843 & / NPART(JEMISS) + UUMNR )**
844 & ( 1.D+00 / (NEXTOT-1) )
845 3800 CONTINUE
846 IF ( LSPREJ ) THEN
847 CALL GRNDM(RNDM,1)
848 IF ( RNDM (1) .GE. FREJE ) GO TO 4100
849 END IF
850 EPS = TVEFF - UURND - BNDGAV
851 4000 CONTINUE
852 IF ( LPRFIX ) THEN
853 ITACC = ITACC + 1
854 EPSFIX (ITACC) = EPS - EPSMIN + BNDGEN
855 END IF
856 DEPSEX = EPS - EPSEEX (JEMISS)
857 IF ( DEPSEX .GT. 0.D+00 ) THEN
858 IF ( EEXDLT (JEMISS) .LT. EEXMNM (JEMISS) .AND.
859 & DEPSEX .GT. 0.5D+00 * ( EPSDLT (JEMISS) - EPSEEX
860 & (JEMISS) ) ) THEN
861 DEPSEX = EPS - EPSDLT (JEMISS)
862 IF ( DEPSEX .GT. 0.5D+00 * ( EPSDLT (JEMISS)
863 & - EPSANY (JEMISS) ) ) THEN
864 EPS = EPSANY (JEMISS)
865 ELSE
866 EPS = EPSDLT (JEMISS)
867 END IF
868 ELSE
869 EPS = EPSEEX (JEMISS)
870 END IF
871 IF ( EPS .LE. EPSMIN ) GO TO 4100
872 END IF
873 AMNHLP = UMO - EPS + EPSMIN - AMNUCL (IEMISS)
874 ERNCM (JEMISS) = 0.5D+00 * ( UMO2 + AMNHLP**2
875 & - AMNUSQ (IEMISS) ) / UMO
876 EPS = UMO - ERNCM (JEMISS) - AMNUCL (IEMISS)
877 EPSINV = 0.5D+00 * ( UMO2 - ( AMNUCL (IEMISS)
878 & + AMNHLP )**2 ) / AMNHLP
879 EPSWLL = EPSINV + AEFRAV + BNDUSE + AMNUCL (IEMISS)
880 BETWLL = SQRT ( 1.D+00 - AMNUSQ (IEMISS) / EPSWLL**2 )
881 EKEWLL = EPSWLL - AMNUCL (IEMISS)
882 EPSCOU = UMO - AMNHLP - AMNUCL (IEMISS)
883 IF ( JEMISS .EQ. 1 ) THEN
884 GSNGLV = 1.5D+00 * ( ANOW - ZNOW ) * SQRT ( ( EKEWLL
885 & + EKFCEN (IEMISS) - AEFRAV ) / EKFCEN (IEMISS))
886 & / EKFCEN (IEMISS)
887 AALPHA = 0.76D+00 + 2.2D+00 / RMASS (IBRES-1)
888 BBETAA = ( 2.22D+00 / RMASS (IBRES-1)**2 - 0.05D+00 )
889 & / AALPHA * 1.D-03
890 SIGINV = 0.1D+00 * XINNEU ( EPSINV, ZAFTR (JEMISS),
891 & BBETAA )
892 IF ( EKEWLL .LE. 0.04D+00 ) THEN
893 SIGMAP = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
894 & -1.86D+00 + 0.09415D+03 * EKEWLL
895 & + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
896 & * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
897 & + 0.13D+03 * EKEWLL )**2 )
898 IF ( EKEWLL .LT. 0.02D+00 ) THEN
899 SIGMAN = 0.3333333333333333D+00 * SIGMAP
900 ELSE
901 SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00
902 & / BETWLL + 42.9D+00
903 END IF
904 ELSE
905 SIGMAP = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
906 & + 82.2D+00
907 SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
908 & + 42.9D+00
909 END IF
910 SIGMAP = 0.1D+00 * SIGMAP
911 SIGMAN = 0.1D+00 * SIGMAN
912 SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
913 & + ( ZNOW - ANPROT ) * SIGMAP )
914 & / ( ANOW - ANPTOT )
915 ELSE
916 GSNGLV = 1.5D+00 * ZNOW * SQRT ( ( EKEWLL
917 & + EKFCEN (IEMISS) - AEFRAV ) / EKFCEN (IEMISS))
918 & / EKFCEN (IEMISS)
919 SIGINV = 1.D-01 * XINPRO ( EPSINV, ZAFTR (JEMISS),
920 & ETHRES )
921 IF ( EKEWLL .LE. 0.04D+00 ) THEN
922 SIGMAN = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
923 & -1.86D+00 + 0.09415D+03 * EKEWLL
924 & + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
925 & * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
926 & + 0.13D+03 * EKEWLL )**2 )
927 IF ( EKEWLL .LT. 0.02D+00 ) THEN
928 SIGMAP = 0.3333333333333333D+00 * SIGMAN
929 ELSE
930 SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00
931 & / BETWLL + 42.9D+00
932 END IF
933 ELSE
934 SIGMAN = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
935 & + 82.2D+00
936 SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
937 & + 42.9D+00
938 END IF
939 SIGMAP = 0.1D+00 * SIGMAP
940 SIGMAN = 0.1D+00 * SIGMAN
941 SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
942 & + ( ZNOW - ANPROT ) * SIGMAP )
943 & / ( ANOW - ANPTOT )
944 END IF
945 ZITA = AEFRAV / EKEWLL
946 IF ( ZITA .LE. 0.5D+00 ) THEN
947 PZITA = 1.D+00 - 1.4D+00 * ZITA
948 ELSE
949 PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
950 & * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
951 END IF
952 RATEC = SIGINV * RATEC0 * EPSINV / GSNGLV
953 IF ( RATEC .LE. 0.D+00 ) GO TO 4100
954 ALAMDC = BETWLL / RATEC
955 ALAMDI = 1.D+00 / ( SGNUNU * RHONUC * PZITA ) + MAX (
956 & TWOTWO * PI * PLABRC * BETWLL * EPSWLL
957 & / AMNUSQ (IEMISS), RZNUCL )
958 AMUTOT = 1.D+00 / ALAMDI + 1.D+00 / ALAMDC
959 RATEI = BETWLL / ALAMDI
960 FLAMDA = RATEC / ( RATEC + CLAMDI * RATEI ) / FLAMMX
961 CALL GRNDM(RNDM,1)
962 IF ( RNDM (1) .LT. FLAMDA ) GO TO 4200
963 4100 CONTINUE
964 IF ( JEMISS .EQ. JEMIS1 ) GO TO 4500
965 IF ( ITACC .LE. NPART (JEMISS) ) THEN
966 NPRFIX (JEMISS) = NPRFIX (JEMISS) + ITACC
967 DO 4150 IT = 1, ITACC
968 EPRFIX = EPRFIX + EPSFIX (IT)
969 4150 CONTINUE
970 ELSE
971 NPRFIX (JEMISS) = NPART (JEMISS)
972 ACCEP = 0.D+00
973 DO 4160 IT = 1, ITACC
974 PRACC = ( NPART (JEMISS) - ACCEP )
975 & / ( ITACC - IT + 1 )
976 CALL GRNDM(RNDM,1)
977 IF ( RNDM (1) .LT. PRACC ) THEN
978 EPRFIX = EPRFIX + EPSFIX (IT)
979 ACCEP = ACCEP + 1.D+00
980 END IF
981 4160 CONTINUE
982 END IF
983 GO TO 4500
984 4200 CONTINUE
985 ZNOW = ZAFTR (JEMISS)
986 ANOW = AAFTR
987 IBRES = IBRES - 1
988 ICRES = ICRES - JEMISS + 1
989 AMMRES = AMMAFT (JEMISS)
990 AMNRES = AMNAFT (JEMISS)
991 EEPCM = EPS + AMNUCL (IEMISS)
992 ERESCM = UMO - EEPCM
993 IF ( LNUCTS ) THEN
994 ELSE IF ( ICYCL .EQ. 1 .AND. LEMISS ) THEN
995 C (1) = C2ND (1)
996 C (2) = C2ND (2)
997 C (3) = C2ND (3)
998 ELSE IF ( ICYCL .LE. 1 .AND. PCMORI .LE. ANGLGB ) THEN
999 C (1) = C1ST (1)
1000 C (2) = C1ST (2)
1001 C (3) = C1ST (3)
1002 ELSE IF ( ICYCL .LE. 1 ) THEN
1003 PCMSXT = - AMNUCL (IEMISS) * ETAX
1004 PCMSYT = - AMNUCL (IEMISS) * ETAY
1005 PCMSZT = - AMNUCL (IEMISS) * ETAZ
1006 APSHAR = MAX ( ANPTOT + NEMISS - 1, ONEONE )
1007 PCMSX = ( CXAXCM * PCMORI + PCMSXT ) / APSHAR - PCMSXT
1008 PCMSY = ( CYAXCM * PCMORI + PCMSYT ) / APSHAR - PCMSYT
1009 PCMSZ = ( CZAXCM * PCMORI + PCMSZT ) / APSHAR - PCMSZT
1010*
1011 PCMSTT = PCMSXT**2 + PCMSYT**2 + PCMSZT**2
1012 PCMSPT = PCMSX**2 + PCMSY**2 + PCMSZ**2
1013 ECMSTT = SQRT ( AMNUSQ (IEMISS) + PCMSTT )
1014 ECMSPT = SQRT ( AMNUSQ (IEMISS) + PCMSPT )
1015*
1016 UMMO2 = 2.D+00 * AMNUSQ (IEMISS) + 2.D+00 * ECMSTT
1017 & * ECMSPT - 2.D+00 * ( PCMSXT * PCMSX + PCMSYT
1018 & * PCMSY + PCMSZT * PCMSZ )
1019 UMMO = SQRT ( UMMO2 )
1020 GAMCCM = ( ECMSTT + ECMSPT ) / UMMO
1021 ETAXM = ( PCMSX + PCMSXT ) / UMMO
1022 ETAYM = ( PCMSY + PCMSYT ) / UMMO
1023 ETAZM = ( PCMSZ + PCMSZT ) / UMMO
1024 ETAPCM = ETAXM * PCMSX + ETAYM * PCMSY + ETAZM * PCMSZ
1025 PHELP = ECMSPT - ETAPCM / ( GAMCCM + 1.D+00 )
1026 PCCMSX = PCMSX - ETAXM * PHELP
1027 PCCMSY = PCMSY - ETAYM * PHELP
1028 PCCMSZ = PCMSZ - ETAZM * PHELP
1029 PCCMSP = SQRT ( PCCMSX**2 + PCCMSY**2 + PCCMSZ**2 )
1030 CALL RACO (C(1),C(2),C(3))
1031 SCAPRO = PCCMSX * C (1) + PCCMSY * C (2) + PCCMSZ * C (3)
1032 IF ( SCAPRO .LT. 0.D+00 ) THEN
1033 C (1) = - C (1)
1034 C (2) = - C (2)
1035 C (3) = - C (3)
1036 END IF
1037 PCCMSX = C(1) * PCCMSP
1038 PCCMSY = C(2) * PCCMSP
1039 PCCMSZ = C(3) * PCCMSP
1040 ETAPCM = ETAXM * PCCMSX + ETAYM * PCCMSY + ETAZM * PCCMSZ
1041 PHELP = 0.5D+00 * UMMO + ETAPCM / ( GAMCCM + 1.D+00 )
1042 PCMSX = PCCMSX + ETAXM * PHELP
1043 PCMSY = PCCMSY + ETAYM * PHELP
1044 PCMSZ = PCCMSZ + ETAZM * PHELP
1045 PCMSPT = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
1046 C (1) = PCMSX / PCMSPT
1047 C (2) = PCMSY / PCMSPT
1048 C (3) = PCMSZ / PCMSPT
1049 ETAPCM = - ETAPCM
1050 PHELP = 0.5D+00 * UMMO + ETAPCM / ( GAMCCM + 1.D+00 )
1051 PCMSX = - PCCMSX + ETAXM * PHELP
1052 PCMSY = - PCCMSY + ETAYM * PHELP
1053 PCMSZ = - PCCMSZ + ETAZM * PHELP
1054 PCMSPT = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
1055 C2ND (1) = PCMSX / PCMSPT
1056 C2ND (2) = PCMSY / PCMSPT
1057 C2ND (3) = PCMSZ / PCMSPT
1058 ELSE
1059 UMULEG (0) = 1.0D+00
1060 UMULEG (1) = 0.6666666666666667D+00
1061 UMULEG (2) = 0.25D+00
1062 ACOLEG (0) = 0.5D+00
1063 ACOLEG (1) = 1.5D+00 * UMULEG (1)**ICYCL
1064 ACOLEG (2) = 2.5D+00 * UMULEG (2)**ICYCL
1065 COSTHE = COSLEG ( ACOLEG )
1066 SINTHE = SQRT ( ( 1.D+00 - COSTHE ) * ( 1.D+00 + COSTHE )
1067 & )
1068 4300 CONTINUE
1069 CALL GRNDM(RNDM,2)
1070 RPHI1 = 2.D+00 * RNDM (1) - 1.D+00
1071 RPHI2 = 2.D+00 * RNDM (2) - 1.D+00
1072 RPHI12 = RPHI1 * RPHI1
1073 RPHI22 = RPHI2 * RPHI2
1074 RSQ = RPHI12 + RPHI22
1075 IF ( RSQ .GT. 1.D+00 ) GO TO 4300
1076 SINPHI = 2.D+00 * RPHI1 * RPHI2 / RSQ
1077 COSPHI = ( RPHI12 - RPHI22 ) / RSQ
1078 SINT02 = ( 1.D+00 - CZAXCM ) * ( 1.D+00 + CZAXCM )
1079 IF ( SINT02 .LT. ANGLSQ ) THEN
1080 C (1) = COSPHI * SINTHE
1081 C (2) = SINPHI * SINTHE
1082 C (3) = CZAXCM * COSTHE
1083 ELSE
1084 SINTH0 = SQRT ( SINT02 )
1085 UPRIME = SINTHE * COSPHI
1086 VPRIME = SINTHE * SINPHI
1087 COSPH0 = CXAXCM / SINTH0
1088 SINPH0 = CYAXCM / SINTH0
1089 C (1) = UPRIME * COSPH0 * CZAXCM - VPRIME * SINPH0
1090 & + COSTHE * CXAXCM
1091 C (2) = UPRIME * SINPH0 * CZAXCM + VPRIME * COSPH0
1092 & + COSTHE * CYAXCM
1093 C (3) = - UPRIME * SINTH0 + COSTHE * CZAXCM
1094 END IF
1095 END IF
1096 PCMS = SQRT ( ( EEPCM - AMNUCL (IEMISS) ) * ( EEPCM
1097 & + AMNUCL (IEMISS) ) )
1098 UMO2 = ( ERESCM - PCMS ) * ( ERESCM + PCMS )
1099 UMO = SQRT (UMO2)
1100 PCMSX = PCMS * C (1)
1101 PCMSY = PCMS * C (2)
1102 PCMSZ = PCMS * C (3)
1103 ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
1104 NP = NP + 1
1105 KPART (NP) = 1 - 7 * ( JEMISS - 2 )
1106 TKI (NP) = GAMCM * EEPCM + ETAPCM - AMNUCL (IEMISS)
1107 PHELP = ETAPCM / ( GAMCM + 1.D+00 ) + EEPCM
1108 WEI (NP) = WEE
1109 PLBPX = PCMSX + ETAX * PHELP
1110 PLBPY = PCMSY + ETAY * PHELP
1111 PLBPZ = PCMSZ + ETAZ * PHELP
1112 PHELP = SQRT ( PLBPX * PLBPX + PLBPY * PLBPY
1113 & + PLBPZ * PLBPZ )
1114 CXR (NP) = PLBPX / PHELP
1115 CYR (NP) = PLBPY / PHELP
1116 CZR (NP) = PLBPZ / PHELP
1117 PLR (NP) = PHELP
1118 IGREYP = IGREYP + JEMISS - 1
1119 IGREYN = IGREYN + 2 - JEMISS
1120 IF ( JEMISS .EQ. 1 ) THEN
1121 ISTORE = IGREYN
1122 ELSE
1123 ISTORE = IGREYP
1124 END IF
1125 IF ( LNUCTS ) THEN
1126 NPCYCL (ISTORE,JEMISS) = -ICYCL
1127 ELSE
1128 NPCYCL (ISTORE,JEMISS) = ICYCL
1129 END IF
1130 EINTR = EINTR + TKI (NP) + AMNUCL (IEMISS)
1131 PXINTR = PXINTR + CXR (NP) * PLR (NP)
1132 PYINTR = PYINTR + CYR (NP) * PLR (NP)
1133 PZINTR = PZINTR + CZR (NP) * PLR (NP)
1134 IBINTR = IBINTR + IBAR ( KPART (NP) )
1135 ICINTR = ICINTR + ICH ( KPART (NP) )
1136 ERES = GAMCM * ERESCM - ETAPCM
1137 EKRES = ERES - UMO
1138 PHELP = - ETAPCM / ( GAMCM + 1.D+00 ) + ERESCM
1139 PXRES = - PCMSX + ETAX * PHELP
1140 PYRES = - PCMSY + ETAY * PHELP
1141 PZRES = - PCMSZ + ETAZ * PHELP
1142 PTRES2= PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
1143 PTRES = SQRT (PTRES2)
1144 TVCMS = UMO - AMNRES
1145 IF ( TVCMS .LT. 0.D+00 ) THEN
1146 IF ( TVCMS .LT. - GAMMIN ) THEN
1147 WRITE (LUNOUT,*)' PREPRE: TVCMS',TVCMS
1148 WRITE (LUNERR,*)' PREPRE: TVCMS',TVCMS
1149 WRITE (77,*)' ^^^ PREPRE: TVCMS',TVCMS
1150 END IF
1151 TVCMS = 0.D+00
1152 END IF
1153 EHLFIX = MIN ( TVCMS, EHLFIX )
1154 TVEFF = TVCMS - EHLFIX
1155 IF ( TVEFF .LE. 0.D+00 ) GO TO 7000
1156 GAMCM = ERES / UMO
1157 ETACM = PTRES / UMO
1158 ETAX = PXRES / UMO
1159 ETAY = PYRES / UMO
1160 ETAZ = PZRES / UMO
1161 ETAPCM = ETAX * PXORI + ETAY * PYORI + ETAZ * PZORI
1162 PHELP = EKORI + AM (KPORI) - ETAPCM / ( GAMCM + 1.D+00 )
1163 PCMSX = PXORI - ETAX * PHELP
1164 PCMSY = PYORI - ETAY * PHELP
1165 PCMSZ = PZORI - ETAZ * PHELP
1166 PCMORI = PCMSX**2 + PCMSY**2 + PCMSZ**2
1167 IF ( PCMORI .LE. ANGLGB ) THEN
1168 PCMORI = 0.D+00
1169 CALL RACO ( CXAXCM, CYAXCM, CZAXCM )
1170 ELSE
1171 PCMORI = SQRT ( PCMORI )
1172 CXAXCM = PCMSX / PCMORI
1173 CYAXCM = PCMSY / PCMORI
1174 CZAXCM = PCMSZ / PCMORI
1175 END IF
1176 IAAFT = IBRES-1
1177 IZAFT = ICRES
1178 AAFTR = IAAFT
1179 ZAFTR (1) = IZAFT
1180 AMMAFT(1) = AAFTR * AMUAMU + 1.D-03
1181 & * FKENER (AAFTR,ZAFTR(1))
1182 AMNAFT(1) = AMMAFT (1) - ZAFTR (1) * AMELEC
1183 & + ELBNDE (IZAFT)
1184 CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (1), EEXMNM (1), EEXDUM )
1185 BNDGAV = BNENRG (2)
1186 EPSMX (1) = UMO - AMNRES - BNDGAV - EHLFIX
1187 EPSDLT (1) = UMO - AMNRES - BNDGAV - EEXDLT (1)
1188 EPSEEX (1) = UMO - AMNRES - BNDGAV - EEXMNM (1)
1189 IZAFT = ICRES - 1
1190 ZAFTR (2) = IZAFT
1191 AMMAFT(2) = AAFTR * AMUAMU + 1.D-03
1192 & * FKENER (AAFTR,ZAFTR(2))
1193 AMNAFT(2) = AMMAFT (2) - ZAFTR (2) * AMELEC
1194 & + ELBNDE (IZAFT)
1195 CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (2), EEXMNM (2), EEXDUM )
1196 BNDGAV = BNENRG (1)
1197 EPSMX (2) = UMO - AMNRES - BNDGAV - EHLFIX
1198 EPSDLT (2) = UMO - AMNRES - BNDGAV - EEXDLT (2)
1199 EPSEEX (2) = UMO - AMNRES - BNDGAV - EEXMNM (2)
1200 EEXANW (1) = 0.D+00
1201 EPSANY (1) = EPSMX (1)
1202 EEXANW (2) = 0.D+00
1203 EPSANY (2) = EPSMX (2)
1204 NPART (JEMISS) = NPART (JEMISS) - 1
1205 IF ( JEMFIX .GT. 0 ) THEN
1206 NPART (1) = NPART (1) + NPRFIX (1)
1207 NPART (2) = NPART (2) + NPRFIX (2)
1208 END IF
1209 NEMISS = NEMISS + 1
1210 LEMISS = .TRUE.
1211 TEMNU = SQRT ( TVCMS / ANOW * ALEVEL )
1212 IF ( EPSMX (1) + EPSMX (2) .LT. 2.D+00 * TEMNU )
1213 & GO TO 7000
1214 SIGIN0 = PI * ( R0SIGK * RMASS (IBRES-1) )**2
1215 NEXMX = SQRT ( 2.D+00 * ANOW * TVCMS / ALEVEL )
1216 ANPROT = NPART (2)
1217 ANNEUT = NPART (1)
1218 NPTOT = NPART (1) + NPART (2)
1219 ANPTOT = NPTOT
1220 GO TO 6000
1221 4500 CONTINUE
1222 4600 CONTINUE
1223 IF ( LNUCTS ) THEN
1224 LEMISS = .TRUE.
1225 GO TO 6000
1226 ELSE
1227 LEMISS = .FALSE.
1228 NPART (1) = NPART (1) + NPRFIX (1)
1229 NPART (2) = NPART (2) + NPRFIX (2)
1230 IF ( JEMFIX .GT. 0 ) THEN
1231 TVEFF = TVCMS - EHLFIX
1232 EPSMX (JEMFIX) = EPSMX0
1233 ERNCM (JEMFIX) = ERNCM0
1234 END IF
1235 END IF
1236 5000 CONTINUE
1237 ANPROT = NPART (2)
1238 ANNEUT = NPART (1)
1239 IF ( NPTOT .GT. 0 ) THEN
1240 PNPROT = ( ZNOW - ANPROT ) * ( 3.D+00 * ANNEUT + ANPROT )
1241 & / ( ANPROT * ( ZNOW - ANPROT + 3.D+00 * ( ANOW
1242 & - ANNEUT - ZNOW ) ) + ANNEUT * ( 3.D+00 * ( ZNOW
1243 & - ANPROT ) + ANOW - ANNEUT - ZNOW ) )
1244 ELSE
1245 PNPROT = ZNOW / ANOW
1246 END IF
1247 CALL GRNDM(RNDM,1)
1248 IF ( RNDM (1) .LT. PNPROT ) THEN
1249 NPART (2) = NPART (2) + 1
1250 ZCOLL = ZCOLL - 1.D+00
1251 ELSE
1252 NPART (1) = NPART (1) + 1
1253 END IF
1254 NPTOT = NPART (1) + NPART (2)
1255 ANPTOT = NPTOT
1256 NHOLE = NHOLE + 1
1257 ACOLL = ACOLL - 1.D+00
1258 AVEBIN = ( ( BBTAR - ZZTAR - ACOLL + ZCOLL ) * BNENRG (2)
1259 & + ( ZZTAR - ZCOLL ) * BNENRG (1) ) / ( BBTAR - ACOLL )
1260 6000 CONTINUE
1261 IF ( LNUCTS ) THEN
1262 RHOIMP = RHOSAV
1263 EKFIMP = EKFSAV
1264 END IF
1265 IF ( NHOLE + NHLFIX + NPTOT .LT. NEXMX .AND.
1266 & NPTOT .LE. NINT (0.5D+00 * ANOW) ) GO TO 1000
1267 7000 CONTINUE
1268*=== End of subroutine prepre =========================================*
1269 RETURN
1270 END
1271#if defined(CERNLIB_HPUX)
1272$OPTIMIZE ON
1273#endif