source: trunk/MagicSoft/Simulation/Corsika/Mmcs/venlnk.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: 16.9 KB
Line 
1 SUBROUTINE VENLNK
2
3C-----------------------------------------------------------------------
4C VEN(US) L(I)NK (TO CORSIKA)
5C
6C LINKS VENUS PACKAGE TO CORSIKA, NEEDS FIRST CALL OF VENINI
7C THIS SUBROUTINE IS CALLED FROM SDPM
8C
9C DESIGN : D. HECK IK3 FZK KARLSRUHE
10C-----------------------------------------------------------------------
11
12*KEEP,INTER.
13 COMMON /INTER/ AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
14 * FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
15 * RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
16 * WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN,
17 * IDIF,ITAR
18 DOUBLE PRECISION AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
19 * FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
20 * RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
21 * WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN
22 INTEGER IDIF,ITAR
23*KEEP,PAM.
24 COMMON /PAM/ PAMA,SIGNUM
25 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
26*KEEP,PARPAR.
27 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
28 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
29 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
30 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
31 INTEGER ITYPE,LEVL
32*KEEP,PARPAE.
33 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
34 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
35 * (CURPAR(4), PHI ), (CURPAR(5), H ),
36 * (CURPAR(6), T ), (CURPAR(7), X ),
37 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
38 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
39 * (CURPAR(12),ECM )
40*KEEP,RANDPA.
41 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
42 DOUBLE PRECISION FAC,U1,U2
43 REAL RD(3000)
44 INTEGER ISEED(103,10),NSEQ
45 LOGICAL KNOR
46*KEEP,REST.
47 COMMON /REST/ CONTNE,TAR,LT
48 DOUBLE PRECISION CONTNE(3),TAR
49 INTEGER LT
50*KEEP,RUNPAR.
51 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
52 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
53 * MONIOU,MDEBUG,NUCNUC,
54 * CETAPE,
55 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
56 * N1STTR,MDBASE,
57 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
58 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
59 * ,GHEISH,GHESIG
60 COMMON /RUNPAC/ DSN,HOST,USER
61 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
62 REAL STEPFC
63 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
64 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
65 * N1STTR,MDBASE
66 INTEGER CETAPE
67 CHARACTER*79 DSN
68 CHARACTER*20 HOST,USER
69
70 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
71 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
72 * ,GHEISH,GHESIG
73*KEEP,VENUS.
74 COMMON /VENUS/ ISH0,IVERVN,MTAR99,FVENUS,FVENSG
75 INTEGER ISH0,IVERVN,MTAR99
76 LOGICAL FVENUS,FVENSG
77*KEND.
78
79 PARAMETER (KOLLMX=2500)
80 PARAMETER (MXEPS=10)
81 PARAMETER (NDEP=129)
82 PARAMETER (NDET=129)
83 PARAMETER (NPRBMS=20)
84 PARAMETER (NPTQ=129)
85 PARAMETER (NSTRU=2049)
86 COMMON /ACCUM/ AMSAC,ILAMAS,IMSG,INOIAC,IPAGE,JERR,NAEVT,NREVT
87 * ,NRPTL,NRSTR,NTEVT
88 COMMON /CDEN/ MASSNR,RMX,R0
89 COMMON /CIPIO/ IPIO
90 COMMON /CNSTA/ AINFIN,PI,PIOM,PROM
91 COMMON /COL/ BIMP,BMAX,COORD(4,KOLLMX),DISTCE(KOLLMX)
92 * ,QDEP(NDEP),QDET14(NDET),QDET16(NDET),QDET40(NDET)
93 * ,QDET99(NDET),RMPROJ,RMTARG(4),XDEP(NDEP)
94 * ,XDET14(NDET),XDET16(NDET),XDET40(NDET)
95 * ,XDET99(NDET)
96 * ,KOLL,LTARG,NORD(KOLLMX),NPROJ,NRPROJ(KOLLMX)
97 * ,NRTARG(KOLLMX),NTARG
98 COMMON /CPRBMS/ PRBMS(NPRBMS)
99 COMMON /CPTQ/ QPTH(NPTQ),QPTQ(NPTQ),XPTQ(NPTQ),QPTQMX,QPTHMX
100 DOUBLE PRECISION SEEDC,SEEDI
101 COMMON /CSEED/ SEEDC,SEEDI
102 COMMON /FILES/ IFCH,IFDT,IFHI,IFMT,IFOP
103 COMMON /NEVNT/ NEVNT
104 COMMON /PARO1/ AMPRIF,AMSIAC,BMAXIM,BMINIM,CORE,CUTMSQ,CUTMSS
105 * ,DELMSS,DELREM,FCTRMX,GAUMX,OVERLP,PAREA,PDIQUA
106 * ,PHARD,PSPINL,PSPINH,PISPN,PTF,PTH,PTMX,PTQ,PUD
107 * ,PVALEN,QSEPC,QSETC,QMUST,QVAPC,QVATC,RADIAC
108 * ,RADIAS,RSTRAS,SIGJ,SIGPPI,TAUMAX,TAUMIN
109 * ,TAUMX,TAUNLL,TENSN,THEMAS,WPROJ,WTARG,WTMINI
110 * ,WTSTEP,XCUT
111 * ,IAQU,IFRADE,IOJINT,IOPBRK,IOPENT,IOPENU
112 * ,IOPTF,IOPTQ,IRESCL,IWCENT,KENTRO,KO1KO2
113 * ,LABSYS,MAXRES,NCLEAN,NCOLMX,NDECAW,NEQMN,NEQMX
114 * ,NSTTAU,NTRYMX,NUMTAU
115 COMMON /PARO2/ AMPROJ,AMTARG,ANGMUE,ELEPTI,ELEPTO,ENGY
116 * ,PNLL,PNLLX,PROB(99),PROSEA,RHOPHI,TAUREA
117 * ,YHAHA,YMXIMI,YPJTL
118 * ,ICBAC(99,2),ICFOR(99,2),ICHOIC,ICLHIS,IDPM
119 * ,IDPROJ,IDTARG,IENTRO,IJPHIS,IMIHIS,IPAGI,ISH
120 * ,ISHEVT,ISHSUB,ISPALL,ISPHIS,ISTMAX,ISUP,IVI
121 * ,JPSI,JPSIFI,KUTDIQ,LAPROJ,LATARG,MAPROJ,MATARG
122 * ,MODSHO,NDECAX,NDECAY,NEVENT
123 COMMON /PARO3/ ASUHAX(7),ASUHAY(7),OMEGA,SIGPPD,SIGPPE,UENTRO
124 * ,IWZZZZ
125 COMMON /PARO4/ GRICEL,GRIDEL,GRIGAM,GRIRSQ,GRISLO
126 COMMON /PARO5/ DELEPS,DELVOL
127 COMMON /QUARKM/ SMAS,SSMAS,USMAS,UUMAS
128 COMMON /STRU/ QSEP(NSTRU),QSET(NSTRU),QVAP(NSTRU)
129 * ,QVAT(NSTRU),XCUTAR,XSTRU(NSTRU)
130 * ,IDTG
131 COMMON /STRU2/ DELTA0,DELTA1,QSEH(NSTRU),QSEPI(NSTRU)
132 * ,QVAH(NSTRU),QVAPI(NSTRU),XSE(NSTRU),XVA(NSTRU)
133
134 DOUBLE PRECISION ERRER,VALUE
135 INTEGER IFLAG
136
137 COMMON /VENLIN/ PTQ1,PTQ2,PTQ3,QMUST1,QMUST2,QMUST3
138 * ,IDTABL(100)
139
140 EXTERNAL SDENSI,SPTQ,SSE0,SVA0,SVA1
141C-----------------------------------------------------------------------
142
143 IF ( DEBUG ) WRITE(MDEBUG,*) 'VENLNK: TAR',SNGL(TAR)
144
145 NSTRUC = NSTRU
146 IF ( DEBUG ) THEN
147 ISH = ISH0
148 ELSE
149 ISH = 0
150 ENDIF
151 NEVNT = SHOWNO
152C SET RANDOM NUMBER GENERATOR STATUS
153 SEEDC=ISEED(2,1)+1.D9*ISEED(3,1)
154C CALCULATE ENERGY IN LAB SYSTEM FOR ELASTICITY FOR VARIOUS PROJECTILES
155 IF ( ITYPE .EQ. 1 ) THEN
156C TREAT PHOTON PROJECTILES (FROM EGS)
157 CALL RMMAR(RD,1,1)
158 IF ( RD(1) .LE. 0.5 ) THEN
159 ITYPE = 7
160 ELSE
161 ITYPE = 17
162 ENDIF
163 ELAB = CURPAR(2)
164 CURPAR(2) = ELAB / PAMA(ITYPE)
165 ELSEIF ( ITYPE .LT. 100 ) THEN
166C TREAT ORDINARY PROJECTILES
167 ELAB = CURPAR(2) * PAMA(ITYPE)
168 ELSE
169C TREAT NUCLEI PROJECTILES
170 NPROT = MOD(ITYPE,100)
171 NNEUT = ITYPE/100 - NPROT
172 ELAB = CURPAR(2) * ( PAMA(14)*NPROT + PAMA(13)*NNEUT )
173 ENDIF
174C SET TARGET PARAMETERS
175 MATARG = NINT(TAR)
176 IDTARG = 1120
177 AMTARG = PAMA(14)
178 IF ( TAR. EQ. 14.D0 ) THEN
179 LTARG = 1
180 LATARG = 7
181 ELSEIF ( TAR .EQ. 16.D0 ) THEN
182 LTARG = 2
183 LATARG = 8
184 ELSEIF ( TAR .EQ. 40.D0 ) THEN
185 LTARG = 3
186 LATARG = 18
187 ELSE
188 WRITE(MONIOU,*)'VENLNK: UNDEFINED TARGET TAR=',SNGL(TAR)
189 ENDIF
190
191C FOR THE CASE OF AN ARBITRARY TARGET (NOT AIR)
192 IF ( LTARG .GT. 3 ) THEN
193 MASSNR = MATARG
194 IF ( MASSNR .GT. 1 ) THEN
195 IF ( MASSNR .NE. MTAR99 ) THEN
196 R0 = 1.19*MASSNR**(.3333333) -1.61*MASSNR**(-.3333333)
197 CX = R0+FCTRMX*0.54
198 RMTARG(4) = CX
199 CALL UTQUAF(SDENSI,NDET,XDET99,QDET99,0.,.33*CX,.66*CX,CX)
200 MTAR99 = MATARG
201 ENDIF
202 ELSE
203 RMTARG(4) = 0.
204 ENDIF
205 ENDIF
206
207C SET PROJECTILE PARAMETERS
208 IF ( ITYPE .LT. 100 ) THEN
209 IDPROJ = IDTABL(ITYPE)
210 IF ( IDPROJ .EQ. 20 .OR. IDPROJ .EQ. -20 ) THEN
211C TREAT NEUTRAL KAONS (K(0)S AND K(0)L)
212 CALL RMMAR(RD,1,1)
213 IF ( RD(1) .LE. 0.5 ) THEN
214 IDPROJ = 230
215 ELSE
216 IDPROJ = -230
217 ENDIF
218 ELSEIF ( IDPROJ .EQ. 2130 ) THEN
219C VENUS CANNOT TREAT LAMBDA, TAKE INSTEAD SIGMA(0))
220 IDPROJ = 1230
221 ELSEIF ( IDPROJ .EQ. -2130 ) THEN
222C VENUS CANNOT TREAT ANTI-LAMBDA, TAKE INSTEAD ANTI-SIGMA(0))
223 IDPROJ = -1230
224 ENDIF
225C ALL OTHER PARTICLE CODES UNCHANGED
226 CALL IDMASS(IDPROJ,AMPROJ)
227 LAPROJ = -1
228 MAPROJ = 1
229 PNLL = CURPAR(2)*AMPROJ
230 ELSE
231C PROJECTILE IS NUCLEUS
232 IDPROJ = 1120
233 CALL IDMASS(IDPROJ,AMPROJ)
234 LAPROJ = MOD(ITYPE,100)
235 MAPROJ = ITYPE/100
236 PNLL = CURPAR(2)*(PAMA(14)+PAMA(13))*0.5
237 ENDIF
238
239 IF ( ABS(IDPROJ) .LT. 1000 ) THEN
240 IF ( ABS(IDPROJ) .EQ. 230 .OR. ABS(IDPROJ) .EQ. 130 ) THEN
241C DIFFRACTIVE PROBABILITY FOR KAON PROJECTILES
242 WPROJ = 0.24
243 ELSE
244C DIFFRACTIVE PROBABILITY FOR PION PROJECTILES
245 WPROJ = 0.20
246 ENDIF
247 ELSE
248C DIFFRACTIVE PROBABILITY FOR BARYON PROJECTILES
249 WPROJ = 0.32
250 ENDIF
251C DIFFRACTIVE PROBABILITY FOR TARGET (ALWAYS NUCLEONS)
252 WTARG = 0.32
253
254 ENGY = SQRT( 2.*SQRT(PNLL**2+AMPROJ**2)*AMTARG+AMTARG**2
255 * +AMPROJ**2 )
256 IF ( DEBUG ) WRITE(MDEBUG,*)'VENLNK: ELAB = ',PNLL,
257 * ' ENGY = ',ENGY
258CDH IF ( ENGY .LT. 12. ) THEN
259 IF ( ENGY .LT. 9.5 ) THEN
260 WRITE (IFMT,*)'VENLNK: ENGY, IDPROJ=',ENGY,IDPROJ
261 CALL UTSTOP('VENLNK: INCIDENT ENERGY TOO SMALL ')
262 ENDIF
263 ENGYI = ENGY
264 PNLLI = PNLL
265 IF ( PNLL .LT. 1.E2 * AMPROJ ) THEN
266 TRM = SQRT(PNLL**2+AMPROJ**2)
267 ENGY = SQRT((TRM+AMTARG-PNLL)*(TRM+AMTARG+PNLL))
268 ELSE
269 TRM = AMPROJ**2*0.5/PNLL+AMTARG
270 ENGY = SQRT(TRM*(2.*PNLL+TRM))
271 ENDIF
272 D1 = ABS(PNLLI-PNLL)/PNLL
273 D2 = ABS(ENGYI-ENGY)/ENGY
274 IF ( D1 .GT. 1.E-3 .OR. D2 .GT. 1.E-3 ) THEN
275 IF ( ISH .GE. 0 ) THEN
276 CALL UTMSG('VENLNK')
277 WRITE(IFCH,*)'***** PNLL,PNLLI:',PNLL,PNLLI
278 WRITE(IFCH,*)'***** ENGY,ENGYI:',ENGY,ENGYI
279 CALL UTMSGF
280 ENDIF
281 ENDIF
282 S = ENGY**2
283 SROOTI = 1./ENGY
284 PNLLX = UTPCM(ENGY,AMPROJ,AMTARG)
285 YHAHA = LOG((SQRT(PNLL**2+S)+PNLL)/ENGY)
286 YPJTL = LOG((SQRT(PNLL**2+AMPROJ**2)+PNLL)/AMPROJ)
287 IF ( ISH .GE. 91 ) WRITE(IFCH,*)'VENLNK: YPJTL=',YPJTL
288
289 ENGYLG = LOG(ENGY)
290 QMUST = QMUST1+QMUST2*ENGYLG+QMUST3*ENGYLG**2
291 PTQ = PTQ1+PTQ2*ENGYLG+PTQ3*ENGYLG**2
292CDH PHARD = 0.030+0.12*(LOG10(S)-LOG10(30.**2))
293 PHARD = 0.030+0.12*(LOG10(S)-2.9542425)
294 PHARD = MIN(1.,PHARD)
295 PHARD = MAX(0.030,PHARD)
296
297C PROJECTILE
298 XCUT = CUTMSQ*SROOTI
299 XCUT2 = XCUT**2
300 IF ( ABS(IDPROJ) .GE. 1000 ) THEN
301C STRUCTURE FUNCTION INTEGRAL FOR BARYONS OF PROJECTILE
302 IPIO = 0
303 CALL UINTEG(VALUE,SSE0,0.D0,1.D0,0.D0,1.D-5,1,ERRER,IFLAG)
304 IF ( IFLAG .GT. 3 .AND. ISH .GT. 0 )
305 * WRITE(IFCH,*)'VENLNK: SSE0:IFLAG=',IFLAG
306 QSEPC = VALUE
307 CALL UINTEG(VALUE,SVA0,0.D0,1.D0,0.D0,1.D-5,1,ERRER,IFLAG)
308 IF ( IFLAG .GT. 3 .AND. ISH .GT. 0 )
309 * WRITE(IFCH,*)'VENLNK: SVA0:IFLAG=',IFLAG
310 QVAPC = VALUE
311 ELSE
312C STRUCTURE FUNCTION INTEGRAL FOR MESONS OF PROJECTILE
313 IPIO = 1
314 A0 = -5.0 + 6.6666667*XCUT2 - 0.53333333*XCUT2**2
315 A1 = 5.0 - 1.875*XCUT2
316 A2 = -3.3333333 + 0.26666667*XCUT2
317 A3 = 1.25
318 A4 = -0.2
319 ROOT = SQRT(XCUT2+1.)
320 QSEPC = 0.9*( (1.-XCUT2*A1)*( LOG(1.+ROOT)-LOG(XCUT) )
321 * - XCUT*A0 + ROOT*(A0+A1+A2+A3+A4) )
322 CALL UINTEG(VALUE,SVA1,0.D0,1.D0,0.D0,1.D-5,1,ERRER,IFLAG)
323 IF ( IFLAG .GT. 3 .AND. ISH .GT. 0 )
324 * WRITE(IFCH,*)'VENLNK: SVA1:IFLAG=',IFLAG
325 QVAPC = VALUE
326 ENDIF
327 IDTG = IPIO
328
329C TARGET
330 IF ( IDTG .EQ. 1 ) THEN
331 IF ( ABS(IDTARG) .GE. 1000 ) THEN
332C STRUCTURE FUNCTION INTEGRAL FOR BARYONS OF TARGET
333 IPIO = 0
334 CALL UINTEG(VALUE,SSE0,0.D0,1.D0,0.D0,1.D-5,1,ERRER,IFLAG)
335 IF ( IFLAG .GT. 3 .AND. ISH .GT. 0 )
336 * WRITE(IFCH,*)'VENLNK: SSE0:IFLAG=',IFLAG
337 QSETC = VALUE
338 CALL UINTEG(VALUE,SVA0,0.D0,1.D0,0.D0,1.D-5,1,ERRER,IFLAG)
339 IF ( IFLAG .GT. 3 .AND. ISH .GT. 0 )
340 * WRITE(IFCH,*)'VENLNK: SVA0:IFLAG=',IFLAG
341 QVATC = VALUE
342 ELSE
343 IPIO=1
344 QVATC = QVAPC
345 QSETC = QSEPC
346 ENDIF
347 ELSE
348 IF ( ABS(IDTARG) .GE. 1000 ) THEN
349 IPIO = 0
350 QVATC = QVAPC
351 QSETC = QSEPC
352 ELSE
353C STRUCTURE FUNCTION INTEGRAL FOR BARYONS OF TARGET
354 IPIO=1
355 A0 = -5.0 + 6.6666667*XCUT2 - 0.53333333*XCUT2**2
356 A1 = 5.0 - 1.875*XCUT2
357 A2 = -3.3333333 + 0.26666667*XCUT2
358 A3 = 1.25
359 A4 = -0.2
360 ROOT = SQRT(XCUT2+1.)
361 QSETC = 0.9*( (1.-XCUT2*A1)*( LOG(1.+ROOT)-LOG(XCUT) )
362 * - XCUT*A0 + ROOT*(A0+A1+A2+A3+A4) )
363 CALL UINTEG(VALUE,SVA1,0.D0,1.D0,0.D0,1.D-5,1,ERRER,IFLAG)
364 IF ( IFLAG .GT. 3 .AND. ISH .GT. 0 )
365 * WRITE(IFCH,*)'VENLNK: SVA1:IFLAG=',IFLAG
366 QVATC = VALUE
367 ENDIF
368 ENDIF
369 IF ( ISH .EQ. 16 .OR. DEBUG ) THEN
370 WRITE(IFCH,301) QVAPC, QSEPC, QVATC, QSETC
371 301 FORMAT(' VENLNK: QVAPC, QSEPC, QVATC, QSETC=',4(F10.7,2X))
372 ENDIF
373
374 IF ( PROSEA .GE. 0. ) THEN
375 QVAPC = 1.0
376 QVATC = 1.0
377 QSEPC = PROSEA
378 QSETC = PROSEA
379 ENDIF
380
381 XCUT = CUTMSS*SROOTI
382 XCUTAR = XCUT
383 B = MIN( 0.05, XCUT*500. )
384 A = MIN( 0.2*B, XCUT*100. )
385 PNLLLG = LOG(PNLL)
386 DELTA0 = EXP(-2.791922 - 0.2091742 * PNLLLG)
387 DELTA1 = EXP(-3.885293 - 0.2029558 * PNLLLG)
388 CALL UTQSEA(A,B,1.)
389 IF ( XCUT .LT. 0.04 ) THEN
390 NEND=1.+REAL(NSTRUC)*2./PI*ACOS(1.-2./PI*ACOS(1.-25.*XCUT))
391 ELSE
392 NEND = NSTRUC
393 ENDIF
394
395 IF ( ABS(IDPROJ) .GE. 1000 ) THEN
396 IPIO = 0
397 DO 203 N = 1,NSTRUC
398 QSEP(N) = QSEH(N)
399 203 CONTINUE
400 DO 2031 N = NEND,NSTRUC
401 QVAP(N) = QVAH(N) - DELTA0
402 2031 CONTINUE
403 ELSE
404 IPIO = 1
405 DO 204 N = 1,NSTRUC
406 QSEP(N) = QSEPI(N)
407 204 CONTINUE
408 DO 2041 N = NEND,NSTRUC
409 QVAP(N) = QVAPI(N) - DELTA1
410 2041 CONTINUE
411 ENDIF
412 CALL UTQVAL(QVAP,NEND)
413
414 IF ( IDTG .EQ. 0 ) THEN
415 IF ( ABS(IDTARG) .GE. 1000 ) THEN
416 IPIO = 0
417 DO 205 N=1,NSTRUC
418 QSET(N) = QSEP(N)
419 QVAT(N) = QVAP(N)
420 205 CONTINUE
421 ELSE
422 IPIO = 1
423 DO 209 N = 1,NSTRUC
424 QSET(N) = QSEPI(N)
425 209 CONTINUE
426 DO 2091 N = NEND,NSTRUC
427 QVAT(N) = QVAPI(N) - DELTA1
428 2091 CONTINUE
429 CALL UTQVAL(QVAT,NEND)
430 ENDIF
431
432 ELSE
433 IF ( ABS(IDTARG) .GE. 1000 ) THEN
434 IPIO = 0
435 DO 210 N = 1,NSTRUC
436 QSET(N) = QSEH(N)
437 210 CONTINUE
438 DO 2101 N = NEND,NSTRUC
439 QVAT(N) = QVAH(N) - DELTA0
440 2101 CONTINUE
441 CALL UTQVAL(QVAT,NEND)
442
443 ELSE
444 IPIO = 1
445 DO 216 N=1,NSTRUC
446 QSET(N) = QSEP(N)
447 QVAT(N) = QVAP(N)
448 216 CONTINUE
449 ENDIF
450 ENDIF
451
452 IF ( ISH .EQ. 21 ) THEN
453 CALL UTHSEA
454 CALL UTSTOP(' VENLNK: ')
455 ENDIF
456
457 QPTHMX = 0.5/PTH**2-PTH**2/(2.*(PTH**2+PTMX**2)**2)
458 IF ( IOPTQ .EQ. 2 ) THEN
459 QPTQMX = 1. - EXP(-PI*PTMX**2/(4.*PTQ**2) )
460 ELSEIF ( IOPTQ .EQ. 3 ) THEN
461 QPTQMX = 1. - PTQ**2/(PTQ**2+PTMX**2)
462 ELSE
463 CX = PTMX
464 CALL UTQUAF(SPTQ,NPTQ,XPTQ,QPTQ,0.,.33*CX,.66*CX,CX)
465 ENDIF
466
467 SIGPPI = -1.0
468C CALCULATE ENERGY DEPENDENT CROSS SECTION FOR BARYONS
469 CALL RACPRO('GRI',QMUST,NPRBMS,PRBMS)
470 IF ( ABS(IDPROJ) .LE. 120 .OR. ABS(IDPROJ) .EQ. 220 ) THEN
471C CROSS SECTION FOR PIONS (OR ETA FOR PHOTONS FROM EGS)
472 SIGPPI = SIGPPI * 0.6667
473 ELSEIF ( ABS(IDPROJ) .EQ. 130 .OR. ABS(IDPROJ) .EQ. 230 ) THEN
474C CROSS SECTION FOR KAONS
475 SIGPPI = SIGPPI * 0.5541
476 ENDIF
477
478 MASSNR = MAPROJ
479 RMPROJ = 0.
480 IF ( MASSNR .GT. 1 ) THEN
481 R0 = 1.19*MASSNR**(.3333333) -1.61*MASSNR**(-.3333333)
482 CX = R0+FCTRMX*0.54
483 RMPROJ = CX
484 CALL UTQUAF(SDENSI,NDEP,XDEP,QDEP,0.,.33*CX,.66*CX,CX)
485 ENDIF
486
487 IF ( IDPM .EQ. 1 ) THEN
488 QSEPC = 0.
489 QSETC = 0.
490 ENDIF
491 BMAX = RMPROJ+RMTARG(LTARG)
492
493 IF ( ISH .GE. 91 ) WRITE(IFCH,*)'VENLNK: AVENUS IS NOW CALLED'
494 CALL AVENUS
495
496C NOW BRING PARTICLES TO CORSIKA STACK
497 CALL VSTORE
498
499 IF ( ISH .GE. 91 ) WRITE(IFCH,*)'VENLNK: (EXIT)'
500 RETURN
501 END
Note: See TracBrowser for help on using the repository browser.