]>
Commit | Line | Data |
---|---|---|
9ef1c2d9 | 1 | c----------------------------------------------------------------------- |
2 | subroutine bjinta(ier) | |
3 | c----------------------------------------------------------------------- | |
4 | c fin. state interactions and decays | |
5 | c----------------------------------------------------------------------- | |
6 | include 'epos.inc' | |
7 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
8 | common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctimel/ntc | |
9 | common/col3/ncol,kolpt | |
10 | double precision ttaun,ttau0,rcproj,rctarg | |
11 | common/cttaun/ttaun /cttau0/ttau0 /geom1/rcproj,rctarg | |
12 | logical go | |
13 | ||
14 | call utpri('bjinta',ish,ishini,4) | |
15 | ||
16 | ier=0 | |
17 | ||
18 | if(ncol.eq.0.and.iappl.eq.2)goto1000 | |
19 | if(nevt.ne.1.or.ifrade.eq.0)goto1000 | |
20 | ||
21 | if(iappl.eq.4)then | |
22 | ctp060829 iacn=0 | |
23 | goto5000 | |
24 | endif | |
25 | ||
26 | if(iappl.eq.1)then | |
27 | tauxx=0.7+0.94*max(radnuc(maproj),radnuc(matarg))/(0.5*engy) | |
28 | else | |
29 | tauxx=0. | |
30 | endif | |
31 | tauzz=max(taumin,tauxx) | |
32 | !print*,'====',taumin,tauxx,tauzz | |
33 | ttaus=dble(tauzz) | |
34 | ttau0=dsqrt(rcproj*rctarg) | |
35 | call jtauin ! initialize hyperbola | |
36 | ||
37 | if(iappl.ne.1)goto 5000 | |
38 | ||
39 | ||
40 | c no-secondary-interactions or parton-ladder-fusion | |
41 | c ------------------------------------------------- | |
42 | if(iorsce.eq.0.and.iorsdf.eq.0.and.iorshh.eq.0 | |
43 | & .or.iorsdf.eq.3)then | |
44 | ctp060829 iacn=0 | |
45 | if(iorsdf.eq.3)then | |
46 | c if nptl already very big, clean up useless particles in cptl list. | |
47 | c (do not use it when gakstr() is called (some information lost) | |
48 | if(nclean.gt.0.and.nptl.gt.mxptl/5)then | |
49 | nptli=maproj+matarg+1 | |
50 | do iii=nptli,nptl | |
51 | go=.true. | |
52 | if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false. | |
53 | if(go.and.mod(istptl(iii),10).ne.0)istptl(iii)=99 | |
54 | enddo | |
55 | nptl0=nptl | |
56 | call utclea(nptli,nptl0) | |
57 | endif | |
58 | nptlbpo=nptl | |
59 | call jintpo !parton-ladder-fusion | |
60 | if(ish.ge.2)call alist('parton-ladder-fusion&',nptlbpo+1,nptl) | |
61 | endif | |
62 | goto5000 | |
63 | else | |
64 | stop'bjinta: not supported any more (310305). ' | |
65 | endif | |
66 | ||
67 | 5000 continue | |
68 | ||
69 | nptlbd=nptl | |
70 | ||
71 | call xSpaceTime | |
72 | ||
73 | if(ifrade.eq.0)goto779 !skip decay | |
74 | if(idecay.eq.0)goto779 !skip decay | |
75 | ||
76 | ||
77 | if(ish.ge.2)call alist('final decay&',0,0) | |
78 | if(iappl.eq.7)then | |
79 | nptli=1 | |
80 | else | |
81 | nptli=maproj+matarg+1 | |
82 | endif | |
83 | np1=nptli | |
84 | 41 np2=nptl | |
85 | nptli=np1 | |
86 | ip=np1-1 | |
87 | do while (ip.lt.np2) | |
88 | ip=ip+1 | |
89 | if(istptl(ip).eq.0)then | |
90 | call hdecas(ip,iret) | |
91 | if(iret.eq.1)goto 1001 | |
92 | if(iret.eq.-1)goto 42 | |
93 | c remove useless particles if not enough space | |
94 | if(nclean.gt.0.and.nptl.gt.mxptl/2)then | |
95 | nnnpt=0 | |
96 | do iii=nptli,ip | |
97 | go=.true. | |
98 | if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false. | |
99 | if(go.and.mod(istptl(iii),10).ne.0)then | |
100 | istptl(iii)=99 | |
101 | nnnpt=nnnpt+1 | |
102 | endif | |
103 | enddo | |
104 | if(nnnpt.gt.mxptl-nptl)then | |
105 | nptl0=nptl | |
106 | call utclea(nptli,nptl0) | |
107 | np2=np2-nnnpt | |
108 | ip=ip-nnnpt | |
109 | nptli=ip | |
110 | endif | |
111 | endif | |
112 | endif | |
113 | 42 continue | |
114 | enddo | |
115 | nptli=max(nptli,np1) | |
116 | np1=np2+1 | |
117 | if(np1.le.nptl)then | |
118 | if(ish.ge.2)then | |
119 | if(ish.ge.3)call alist('partial list',0,0) | |
120 | do 6 ip=np1,nptl | |
121 | call alist('&',ip,ip) | |
122 | 6 continue | |
123 | endif | |
124 | goto 41 | |
125 | endif | |
126 | 779 continue | |
127 | ||
128 | if(ish.ge.2)call alist('complete list&',1,nptl) | |
129 | ||
130 | c on shell check | |
131 | c -------------- | |
132 | c if(iappl.eq.1)call jresc | |
133 | ||
134 | 1000 continue | |
135 | call utprix('bjinta',ish,ishini,4) | |
136 | return | |
137 | ||
138 | 1001 continue | |
139 | ier=1 | |
140 | goto1000 | |
141 | ||
142 | end | |
143 | ||
144 | cc---------------------------------------------------------------------- | |
145 | c subroutine jintcs(i,j,ecm,bij,nq,jc,ics) | |
146 | cc---------------------------------------------------------------------- | |
147 | cc compare hadron distance with energy dependent cross section | |
148 | cc data taken from particle data group, durham and juelich | |
149 | cc input: | |
150 | cc i,j: particle indices | |
151 | cc ecm: center-of-mass energy | |
152 | cc bij: impact parameter | |
153 | cc nq: net quark number of fused object | |
154 | cc jc: jc of fused object | |
155 | cc output: | |
156 | cc ics=0 if distance larger than sqrt(sig(E_CMS)/pi) | |
157 | cc ics=1 else | |
158 | cc The data are from HEPDATA, | |
159 | cc the formulas from Rev. Particle Properties 1995 | |
160 | cc---------------------------------------------------------------------- | |
161 | c include 'epos.inc' | |
162 | c integer jci(nflav,2),jcj(nflav,2),jc(nflav,2),kc(nflav) | |
163 | c *,kci(nflav),kcj(nflav) | |
164 | c common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl) | |
165 | c *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3) | |
166 | c common/cdfptl/idfptl(mxptl) | |
167 | c parameter(npp=249,napp=205,npn=411,napn=31,npip=441) | |
168 | c parameter(npim=578,nkmp=299,nkmn=41,nkpp=172,nkpn=91) | |
169 | cc parameter(npim=578,nkmp=299,nkmn=41,nkpp=172,nkpn=91,nlp=35) | |
170 | c parameter(npi1=12,npi2=12,npi3=18,npi4=21,npi5=9) | |
171 | cc real ppecm(npp) | |
172 | c real ppbmx(npp) | |
173 | c real appecm(napp),appbmx(napp) | |
174 | c real pnecm(npn),pnbmx(npn) | |
175 | c real apnecm(napn),apnbmx(napn) | |
176 | c real pipecm(npip),pipbmx(npip) | |
177 | c real pimecm(npim),pimbmx(npim) | |
178 | c real kmpecm(nkmp),kmpbmx(nkmp) | |
179 | c real kmnecm(nkmn),kmnbmx(nkmn) | |
180 | c real kppecm(nkpp),kppbmx(nkpp) | |
181 | c real kpnecm(nkpn),kpnbmx(nkpn) | |
182 | cc real lpecm(nlp),lpbmx(nlp) | |
183 | c real pi1ecm(npi1),pi1bmx(npi1) | |
184 | c real pi2ecm(npi2),pi2bmx(npi2) | |
185 | c real pi3ecm(npi3),pi3bmx(npi3) | |
186 | c real pi4ecm(npi4),pi4bmx(npi4) | |
187 | c real pi5ecm(npi5),pi5bmx(npi5) | |
188 | c | |
189 | cc data ppecm/ | |
190 | cc * 1.8812, 1.8855, 1.8910, 1.8963, 1.9073 | |
191 | cc *, 1.9108, 1.9145, 1.9224, 1.9244, 1.9352 | |
192 | cc *, 1.9466, 1.9468, 1.9542, 1.9592, 1.9636 | |
193 | cc *, 1.9772, 1.9860, 1.9945, 2.0032, 2.0052 | |
194 | cc *, 2.0070, 2.0272, 2.0275, 2.0302, 2.0333 | |
195 | cc *, 2.0402, 2.0427, 2.0586, 2.0608, 2.0692 | |
196 | cc *, 2.0702, 2.0708, 2.0715, 2.0718, 2.0751 | |
197 | cc *, 2.0797, 2.0813, 2.0829, 2.0843, 2.0846 | |
198 | cc *, 2.0935, 2.1062, 2.1113, 2.1123, 2.1170 | |
199 | cc *, 2.1173, 2.1184, 2.1268, 2.1289, 2.1357 | |
200 | cc *, 2.1467, 2.1511, 2.1522, 2.1553, 2.1618 | |
201 | cc *, 2.1639, 2.1726, 2.1771, 2.1795, 2.1799 | |
202 | cc *, 2.1802, 2.1813, 2.1868, 2.2152, 2.2184 | |
203 | cc *, 2.2212, 2.2254, 2.2395, 2.2405, 2.2532 | |
204 | cc *, 2.2606, 2.2613, 2.2861, 2.2889, 2.2914 | |
205 | cc *, 2.2988, 2.3101, 2.3136, 2.3228, 2.3348 | |
206 | cc *, 2.3490, 2.3525, 2.3529, 2.3800, 2.3842 | |
207 | cc *, 2.3912, 2.4088, 2.4130, 2.4193, 2.4298 | |
208 | cc *, 2.4315, 2.4392, 2.4472, 2.4573, 2.5034 | |
209 | cc *, 2.5131, 2.5268, 2.5743, 2.5848, 2.5916 | |
210 | cc *, 2.6327, 2.6620, 2.6700, 2.6984, 2.7080 | |
211 | cc *, 2.7205, 2.7534, 2.7573, 2.7651, 2.7670 | |
212 | cc *, 2.7844, 2.8024, 2.8092, 2.8127, 2.8533 | |
213 | cc *, 2.8556, 2.8638, 2.9079, 2.9395, 2.9500 | |
214 | cc *, 2.9776, 2.9961, 3.0495, 3.0769, 3.0879 | |
215 | cc *, 3.1358, 3.1547, 3.2251, 3.2371, 3.3020 | |
216 | cc *, 3.3527, 3.3620, 3.4221, 3.4967, 3.5019 | |
217 | cc *, 3.5035, 3.5814, 3.5829, 3.6266, 3.8549 | |
218 | cc *, 3.8742, 4.0503, 4.0663, 4.0698, 4.0732 | |
219 | cc *, 4.0778, 4.1074, 4.1300, 4.4976, 4.5183 | |
220 | cc *, 4.5389, 4.5410, 4.5615, 4.6808, 4.9146 | |
221 | cc *, 4.9336, 5.0088, 5.2993, 5.3011, 5.4731 | |
222 | cc *, 5.6083, 5.6416, 5.6465, 5.9171, 5.9502 | |
223 | cc *, 5.9644, 6.1697, 6.1803, 6.2706, 6.3034 | |
224 | cc *, 6.3390, 6.5627, 6.7040, 6.8424, 6.9105 | |
225 | cc *, 6.9780, 7.1111, 7.1662, 7.6202, 7.8624 | |
226 | cc *, 8.2124, 8.7647, 9.0282, 9.2843, 9.5825 | |
227 | cc *, 9.7763, 9.9851, 10.2447, 10.6927, 10.8926 | |
228 | cc *, 11.4549, 11.5365, 11.7779, 11.8519, 13.6241 | |
229 | cc *, 13.6883, 13.7611, 13.8968, 15.0628, 15.1868 | |
230 | cc *, 16.6595, 16.8275, 17.9077, 18.0121, 18.1677 | |
231 | cc *, 19.2213, 19.4156, 19.6556, 19.7002, 21.2604 | |
232 | cc *, 22.9574, 23.3624, 23.4057, 23.4965, 23.4967 | |
233 | cc *, 23.5964, 23.7605, 23.8787, 24.1521, 25.2904 | |
234 | cc *, 26.3796, 27.5960, 30.5240, 30.5954, 30.5957 | |
235 | cc *, 30.6555, 30.6954, 30.7954, 35.1947, 44.6933 | |
236 | cc *, 44.6937, 44.6938, 44.7706, 44.8228, 44.8933 | |
237 | cc *, 45.1933, 52.6798, 52.7090, 52.7921, 52.7921 | |
238 | cc *, 52.7927, 52.8927, 53.1921, 62.2907, 62.3914 | |
239 | cc *, 62.4907, 62.6907, 62.6913, 62.7906 | |
240 | cc */ | |
241 | c data ppbmx/ | |
242 | c * 3.1615, 2.2212, 1.7113, 1.4927, 1.1631 | |
243 | c *, 1.0911, 1.0388, .9525, .9390, .8885 | |
244 | c *, .8019, .8956, .9115, .8593, .8840 | |
245 | c *, .9062, .8444, .8444, .8482, .8713 | |
246 | c *, .8575, .8795, .8795, .8759, .8630 | |
247 | c *, .8822, .8593, .8813, .9036, .9150 | |
248 | c *, .8740, .9219, .8740, .9253, .8722 | |
249 | c *, .8795, .9288, .8704, .9398, .9271 | |
250 | c *, .9407, .9373, .9756, .9926, 1.0357 | |
251 | c *, 1.0037, 1.0408, .9739, 1.0108, 1.0479 | |
252 | c *, 1.0645, 1.0705, 1.1113, 1.0794, 1.0955 | |
253 | c *, 1.1085, 1.1256, 1.1535, 1.1731, 1.1424 | |
254 | c *, 1.0794, 1.1480, 1.1590, 1.1888, 1.2218 | |
255 | c *, 1.2164, 1.2087, 1.2296, 1.2231, 1.2335 | |
256 | c *, 1.2322, 1.2309, 1.2114, 1.2296, 1.2293 | |
257 | c *, 1.2489, 1.2303, 1.2322, 1.2322, 1.2127 | |
258 | c *, 1.2192, 1.2295, 1.2399, 1.2290, 1.2374 | |
259 | c *, 1.2205, 1.2278, 1.2284, 1.2061, 1.2296 | |
260 | c *, 1.2296, 1.2540, 1.2008, 1.2260, 1.2229 | |
261 | c *, 1.2257, 1.2188, 1.2118, 1.2078, 1.1982 | |
262 | c *, 1.2039, 1.2012, 1.1991, 1.1808, 1.1969 | |
263 | c *, 1.1959, 1.1922, 1.1902, 1.1897, 1.1878 | |
264 | c *, 1.1888, 1.1860, 1.1856, 1.1850, 1.2244 | |
265 | c *, 1.1782, 1.1790, 1.1718, 1.1696, 1.1726 | |
266 | c *, 1.1576, 1.1656, 1.1606, 1.1603, 1.1581 | |
267 | c *, 1.1590, 1.1530, 1.1576, 1.1487, 1.1476 | |
268 | c *, 1.1447, 1.1794, 1.1448, 1.1507, 1.1507 | |
269 | c *, 1.1407, 1.1403, 1.1507, 1.1581, 1.1645 | |
270 | c *, 1.1616, 1.1507, 1.1332, 1.1294, 1.1284 | |
271 | c *, 1.1227, 1.1284, 1.1298, 1.1261, 1.1199 | |
272 | c *, 1.1365, 1.1438, 1.1284, 1.1424, 1.1230 | |
273 | c *, 1.1218, 1.1142, 1.1156, 1.1202, 1.1198 | |
274 | c *, 1.1099, 1.1099, 1.1175, 1.1241, 1.1168 | |
275 | c *, 1.1099, 1.1128, 1.1241, 1.1103, 1.1149 | |
276 | c *, 1.1155, 1.1083, 1.1197, 1.1127, 1.1185 | |
277 | c *, 1.1113, 1.1128, 1.1113, 1.1083, 1.1056 | |
278 | c *, 1.1067, 1.1070, 1.1059, 1.1063, 1.1070 | |
279 | c *, 1.1028, 1.0979, 1.1060, 1.1062, 1.1774 | |
280 | c *, 1.0805, 1.1039, 1.0940, 1.0955, 1.0940 | |
281 | c *, 1.0955, 1.1082, 1.1128, 1.1082, 1.0852 | |
282 | c *, 1.1003, 1.1097, 1.1118, 1.0911, 1.1227 | |
283 | c *, 1.1070, 1.1206, 1.1142, 1.1213, 1.1174 | |
284 | c *, 1.1199, 1.1142, 1.1185, 1.1180, 1.1113 | |
285 | c *, 1.1099, 1.1297, 1.1142, 1.1226, 1.1240 | |
286 | c *, 1.1251, 1.1368, 1.1403, 1.1319, 1.1296 | |
287 | c *, 1.1333, 1.1303, 1.1284, 1.1343, 1.1516 | |
288 | c *, 1.1521, 1.1549, 1.1641, 1.1631, 1.1562 | |
289 | c *, 1.1631, 1.1697, 1.1726, 1.1756, 1.1659 | |
290 | c *, 1.1660, 1.1617, 1.1686, 1.1774, 1.1713 | |
291 | c *, 1.1744, 1.1810, 1.1694, 1.1848 | |
292 | c */ | |
293 | c data appecm/ | |
294 | c * 1.9002, 1.9050, 1.9072, 1.9078, 1.9091 | |
295 | c *, 1.9129, 1.9157, 1.9162, 1.9174, 1.9176 | |
296 | c *, 1.9180, 1.9195, 1.9201, 1.9224, 1.9226 | |
297 | c *, 1.9246, 1.9252, 1.9255, 1.9257, 1.9271 | |
298 | c *, 1.9282, 1.9293, 1.9301, 1.9310, 1.9319 | |
299 | c *, 1.9328, 1.9334, 1.9345, 1.9359, 1.9370 | |
300 | c *, 1.9372, 1.9384, 1.9393, 1.9398, 1.9407 | |
301 | c *, 1.9426, 1.9430, 1.9433, 1.9452, 1.9454 | |
302 | c *, 1.9473, 1.9485, 1.9495, 1.9500, 1.9510 | |
303 | c *, 1.9515, 1.9547, 1.9559, 1.9562, 1.9579 | |
304 | c *, 1.9610, 1.9615, 1.9644, 1.9680, 1.9718 | |
305 | c *, 1.9755, 1.9788, 1.9829, 1.9871, 1.9911 | |
306 | c *, 1.9954, 1.9994, 2.0813, 2.0979, 2.1146 | |
307 | c *, 2.1180, 2.1316, 2.1487, 2.1660, 2.1834 | |
308 | c *, 2.1868, 2.1938, 2.1991, 2.2008, 2.2184 | |
309 | c *, 2.2226, 2.2359, 2.2500, 2.2606, 2.2712 | |
310 | c *, 2.2889, 2.2995, 2.3066, 2.3243, 2.3419 | |
311 | c *, 2.3490, 2.3511, 2.3596, 2.3701, 2.3772 | |
312 | c *, 2.3860, 2.3877, 2.3947, 2.4035, 2.4123 | |
313 | c *, 2.4298, 2.4472, 2.4629, 2.4820, 2.4993 | |
314 | c *, 2.5165, 2.5268, 2.5337, 2.5508, 2.5678 | |
315 | c *, 2.5848, 2.5950, 2.6017, 2.6186, 2.6353 | |
316 | c *, 2.6377, 2.6520, 2.6654, 2.6687, 2.6852 | |
317 | c *, 2.7017, 2.7182, 2.7345, 2.7508, 2.7670 | |
318 | c *, 2.7832, 2.7896, 2.7992, 2.8152, 2.8312 | |
319 | c *, 2.8439, 2.8470, 2.8565, 2.8628, 2.9377 | |
320 | c *, 2.9561, 2.9745, 3.0351, 3.0769, 3.0828 | |
321 | c *, 3.1648, 3.2507, 3.2788, 3.3620, 3.4568 | |
322 | c *, 3.5492, 3.6266, 3.7893, 3.8501, 3.8597 | |
323 | c *, 3.8742, 3.9455, 4.1074, 4.3499, 4.3926 | |
324 | c *, 4.5389, 4.5799, 4.6808, 4.7991, 4.9336 | |
325 | c *, 4.9901, 5.1742, 5.2993, 5.3520, 5.4731 | |
326 | c *, 5.6416, 5.6911, 5.8534, 5.9644, 6.0113 | |
327 | c *, 6.2706, 6.3153, 6.7040, 6.9780, 7.3061 | |
328 | c *, 7.6202, 7.7422, 7.8624, 7.8743, 8.0393 | |
329 | c *, 8.2124, 8.4930, 8.7647, 9.0282, 9.2843 | |
330 | c *, 9.5335, 9.7763, 11.5365, 13.7611, 13.8630 | |
331 | c *, 15.0628, 16.8275, 17.9077, 19.4156, 21.2604 | |
332 | c *, 22.9574, 30.4098, 30.5943, 30.6861, 52.5843 | |
333 | c *, 52.7979, 52.7979, 62.2853, 62.4957, 62.6905 | |
334 | c *, 539.9198, 546.9191, 899.8658, 900.0000,1803.0007 | |
335 | c */ | |
336 | c data appbmx/ | |
337 | c * 2.7553, 2.6857, 2.6732, 2.6517, 2.6306 | |
338 | c *, 2.5376, 2.4398, 2.4779, 2.5212, 2.5143 | |
339 | c *, 2.4573, 2.4799, 2.4495, 2.4286, 2.4521 | |
340 | c *, 2.3683, 2.4463, 2.4240, 2.4489, 2.3743 | |
341 | c *, 2.4162, 2.3561, 2.3903, 2.3796, 2.3221 | |
342 | c *, 2.3689, 2.4181, 2.3221, 2.3823, 2.2987 | |
343 | c *, 2.3337, 2.3568, 2.3487, 2.2911, 2.3344 | |
344 | c *, 2.2701, 2.3097, 2.3207, 2.2631, 2.2603 | |
345 | c *, 2.2666, 2.2412, 2.2319, 2.3097, 2.2327 | |
346 | c *, 2.2104, 2.2104, 2.2645, 2.2369, 2.1880 | |
347 | c *, 2.1756, 2.2162, 2.1491, 2.1484, 2.1365 | |
348 | c *, 2.1320, 2.1027, 2.1065, 2.0913, 2.0776 | |
349 | c *, 2.0668, 2.0883, 1.9333, 1.9098, 1.8851 | |
350 | c *, 1.8932, 1.8721, 1.8623, 1.8520, 1.8409 | |
351 | c *, 1.8678, 1.8575, 1.8325, 1.8429, 1.8088 | |
352 | c *, 1.7698, 1.7941, 1.7864, 1.7814, 1.7736 | |
353 | c *, 1.7645, 1.7623, 1.7576, 1.7523, 1.7445 | |
354 | c *, 1.7419, 1.7019, 1.7342, 1.7311, 1.7271 | |
355 | c *, 1.7210, 1.7182, 1.7161, 1.7119, 1.7054 | |
356 | c *, 1.6947, 1.6816, 1.6780, 1.6678, 1.6599 | |
357 | c *, 1.6509, 1.6328, 1.6449, 1.6396, 1.6319 | |
358 | c *, 1.6270, 1.5761, 1.6043, 1.6120, 1.6069 | |
359 | c *, 1.5978, 1.6018, 1.5891, 1.5948, 1.5905 | |
360 | c *, 1.5849, 1.5778, 1.5736, 1.5679, 1.5785 | |
361 | c *, 1.5583, 1.5476, 1.5519, 1.5467, 1.5416 | |
362 | c *, 1.5233, 1.5368, 1.5492, 1.5313, 1.4895 | |
363 | c *, 1.5579, 1.5107, 1.4680, 1.4725, 1.4586 | |
364 | c *, 1.3889, 1.4799, 1.4472, 1.4371, 1.3576 | |
365 | c *, 1.4228, 1.3922, 1.3762, 1.3854, 1.3669 | |
366 | c *, 1.4161, 1.3623, 1.3638, 1.3530, 1.3762 | |
367 | c *, 1.3273, 1.2866, 1.2989, 1.3159, 1.2828 | |
368 | c *, 1.3086, 1.2797, 1.2704, 1.2940, 1.2741 | |
369 | c *, 1.2514, 1.2540, 1.2803, 1.2653, 1.2603 | |
370 | c *, 1.2386, 1.2192, 1.2283, 1.2231, 1.2140 | |
371 | c *, 1.2048, 1.2114, 1.2099, 1.2052, 1.2048 | |
372 | c *, 1.1984, 1.1928, 1.1901, 1.1902, 1.1888 | |
373 | c *, 1.1848, 1.1769, 1.1706, 1.1579, 1.1608 | |
374 | c *, 1.1521, 1.1534, 1.1520, 1.1495, 1.1549 | |
375 | c *, 1.1550, 1.1580, 1.1672, 1.1562, 1.1743 | |
376 | c *, 1.1859, 1.1950, 1.1851, 1.1821, 1.1987 | |
377 | c *, 1.4740, 1.4037, 1.4483, 1.4417, 1.5149 | |
378 | c */ | |
379 | c data pnecm/ | |
380 | c * 1.87867, 1.87869, 1.87872, 1.87875, 1.87878 | |
381 | c *, 1.87881, 1.87884, 1.87887, 1.87890, 1.87893 | |
382 | c *, 1.87896, 1.87899, 1.87902, 1.87906, 1.87909 | |
383 | c *, 1.87913, 1.87916, 1.87920, 1.87923, 1.87927 | |
384 | c *, 1.87931, 1.87934, 1.87938, 1.87942, 1.87946 | |
385 | c *, 1.87950, 1.87954, 1.87958, 1.87962, 1.87966 | |
386 | c *, 1.87970, 1.87975, 1.87979, 1.87983, 1.87988 | |
387 | c *, 1.87992, 1.87997, 1.88001, 1.88006, 1.88011 | |
388 | c *, 1.88015, 1.88020, 1.88025, 1.88030, 1.88035 | |
389 | c *, 1.88040, 1.88045, 1.88050, 1.88055, 1.88061 | |
390 | c *, 1.88066, 1.88071, 1.88077, 1.88082, 1.88087 | |
391 | c *, 1.88093, 1.88099, 1.88104, 1.88110, 1.88116 | |
392 | c *, 1.88121, 1.88127, 1.88133, 1.88139, 1.88145 | |
393 | c *, 1.88151, 1.88157, 1.88163, 1.88170, 1.88176 | |
394 | c *, 1.88182, 1.88189, 1.88195, 1.88202, 1.88208 | |
395 | c *, 1.88215, 1.88221, 1.88228, 1.88235, 1.88241 | |
396 | c *, 1.88248, 1.88255, 1.88262, 1.88269, 1.88276 | |
397 | c *, 1.88283, 1.88290, 1.88297, 1.88305, 1.88312 | |
398 | c *, 1.88319, 1.88327, 1.88334, 1.88342, 1.88349 | |
399 | c *, 1.88357, 1.88364, 1.88372, 1.88380, 1.88388 | |
400 | c *, 1.88396, 1.88403, 1.88411, 1.88419, 1.88428 | |
401 | c *, 1.88436, 1.88444, 1.88452, 1.88460, 1.88469 | |
402 | c *, 1.88477, 1.88485, 1.88494, 1.88502, 1.88511 | |
403 | c *, 1.88519, 1.88528, 1.88537, 1.88546, 1.88554 | |
404 | c *, 1.88563, 1.88572, 1.88581, 1.88590, 1.88599 | |
405 | c *, 1.88608, 1.88618, 1.88627, 1.88636, 1.88645 | |
406 | c *, 1.88655, 1.88664, 1.88674, 1.88683, 1.88693 | |
407 | c *, 1.88702, 1.88712, 1.88722, 1.88731, 1.88741 | |
408 | c *, 1.88751, 1.88761, 1.88771, 1.88781, 1.88791 | |
409 | c *, 1.88801, 1.88811, 1.88822, 1.88832, 1.88842 | |
410 | c *, 1.88852, 1.88863, 1.88873, 1.88884, 1.88894 | |
411 | c *, 1.88905, 1.88916, 1.88926, 1.88937, 1.88948 | |
412 | c *, 1.88959, 1.88970, 1.88980, 1.88991, 1.89003 | |
413 | c *, 1.89014, 1.89025, 1.89036, 1.89047, 1.89058 | |
414 | c *, 1.89070, 1.89081, 1.89093, 1.89104, 1.89116 | |
415 | c *, 1.89127, 1.89139, 1.89150, 1.89162, 1.89174 | |
416 | c *, 1.89186, 1.89198, 1.89209, 1.89221, 1.89233 | |
417 | c *, 1.89245, 1.89258, 1.89270, 1.89282, 1.89294 | |
418 | c *, 1.89306, 1.89319, 1.89331, 1.89344, 1.89356 | |
419 | c *, 1.89369, 1.89381, 1.89394, 1.89406, 1.89419 | |
420 | c *, 1.89432, 1.89445, 1.89458, 1.89496, 1.89549 | |
421 | c *, 1.89602, 1.89615, 1.89656, 1.89697, 1.89724 | |
422 | c *, 1.89766, 1.89780, 1.89850, 1.89893, 1.89922 | |
423 | c *, 1.89995, 1.90068, 1.90098, 1.90143, 1.90174 | |
424 | c *, 1.90235, 1.90250, 1.90312, 1.90407, 1.90503 | |
425 | c *, 1.90519, 1.90616, 1.90732, 1.90749, 1.90833 | |
426 | c *, 1.90936, 1.91216, 1.91379, 1.91545, 1.91714 | |
427 | c *, 1.91905, 1.91924, 1.92100, 1.92160, 1.92179 | |
428 | c *, 1.92259, 1.92319, 1.92421, 1.92502, 1.92543 | |
429 | c *, 1.92605, 1.92646, 1.92792, 1.93046, 1.93241 | |
430 | c *, 1.93306, 1.93570, 1.93704, 1.93953, 1.94021 | |
431 | c *, 1.94183, 1.94699, 1.94723, 1.95206, 1.95329 | |
432 | c *, 1.95452, 1.95651, 1.96029, 1.96080, 1.97927 | |
433 | c *, 1.98533, 1.99806, 2.00884, 2.01269, 2.02779 | |
434 | c *, 2.02963, 2.04329, 2.04643, 2.05915, 2.05947 | |
435 | c *, 2.05979, 2.07207, 2.07337, 2.07533, 2.09178 | |
436 | c *, 2.10847, 2.11385, 2.12061, 2.12536, 2.14243 | |
437 | c *, 2.15308, 2.15343, 2.15964, 2.17072, 2.17698 | |
438 | c *, 2.18185, 2.19442, 2.21193, 2.21369, 2.22353 | |
439 | c *, 2.22846, 2.22952, 2.24715, 2.25597, 2.28851 | |
440 | c *, 2.29134, 2.29382, 2.31293, 2.31328, 2.32673 | |
441 | c *, 2.34935, 2.35501, 2.36207, 2.38251, 2.39729 | |
442 | c *, 2.41135, 2.41591, 2.42817, 2.44006, 2.45995 | |
443 | c *, 2.48220, 2.49952, 2.50643, 2.52951, 2.56750 | |
444 | c *, 2.58756, 2.63546, 2.64719, 2.66487, 2.67285 | |
445 | c *, 2.71088, 2.72336, 2.75178, 2.75634, 2.76802 | |
446 | c *, 2.76867, 2.76996, 2.78741, 2.80542, 2.81567 | |
447 | c *, 2.85639, 2.85860, 2.86681, 2.91102, 2.94265 | |
448 | c *, 2.96543, 3.05273, 3.08729, 3.09115, 3.15805 | |
449 | c *, 3.16821, 3.22857, 3.24053, 3.29551, 3.35628 | |
450 | c *, 3.42580, 3.50727, 3.55297, 3.58520, 3.58675 | |
451 | c *, 3.59580, 3.63049, 3.72425, 3.75633, 4.05460 | |
452 | c *, 4.06219, 4.07412, 4.11175, 4.37426, 4.52314 | |
453 | c *, 4.54378, 4.74537, 4.93885, 5.12513, 5.30495 | |
454 | c *, 5.40999, 5.47892, 5.61424, 5.64757, 5.93444 | |
455 | c *, 5.97071, 6.27732, 6.42517, 6.51227, 6.56970 | |
456 | c *, 6.71114, 6.87703, 6.98545, 7.18446, 7.23904 | |
457 | c *, 7.24942, 7.62830, 8.10603, 8.22114, 8.54948 | |
458 | c *, 8.77406, 9.29418, 9.78673, 10.15712, 10.25566 | |
459 | c *, 10.70407, 11.13446, 11.54882, 12.33587, 13.77577 | |
460 | c *, 15.07880, 15.74960, 16.84544, 17.92675, 18.18703 | |
461 | c *, 18.44365, 19.43624, 20.14863, 21.28302, 22.69375 | |
462 | c *, 22.98187 | |
463 | c */ | |
464 | c data pnbmx/ | |
465 | c * 10.7613, 10.6195, 10.5290, 10.4010, 10.2901 | |
466 | c *, 10.1970, 10.1029, 10.0088, 9.9009, 9.7973 | |
467 | c *, 9.7125, 9.6198, 9.5613, 9.4281, 9.3628 | |
468 | c *, 9.2801, 9.2068, 9.1103, 9.0239, 8.9609 | |
469 | c *, 8.8562, 8.8011, 8.7574, 8.6332, 8.5487 | |
470 | c *, 8.4788, 8.4157, 8.3606, 8.2746, 8.2604 | |
471 | c *, 8.1583, 8.0778, 8.0127, 7.9567, 7.8691 | |
472 | c *, 7.8176, 7.8190, 7.7028, 7.6634, 7.5956 | |
473 | c *, 7.5808, 7.4956, 7.4144, 7.3600, 7.3071 | |
474 | c *, 7.2947, 7.1947, 7.1389, 7.1333, 7.0404 | |
475 | c *, 6.9856, 6.9468, 6.9171, 6.8369, 6.7869 | |
476 | c *, 6.7829, 6.6804, 6.6483, 6.5982, 6.5464 | |
477 | c *, 6.5060, 6.4541, 6.4539, 6.3684, 6.3219 | |
478 | c *, 6.2777, 6.2832, 6.1881, 6.1494, 6.1048 | |
479 | c *, 6.0642, 6.0207, 6.0255, 5.9756, 5.9161 | |
480 | c *, 5.8717, 5.8276, 5.7878, 5.7481, 5.7319 | |
481 | c *, 5.6977, 5.6492, 5.6142, 5.5796, 5.5428 | |
482 | c *, 5.4939, 5.4659, 5.4549, 5.3933, 5.3601 | |
483 | c *, 5.3311, 5.3037, 5.2760, 5.2429, 5.2381 | |
484 | c *, 5.2384, 5.1470, 5.1112, 5.0812, 5.0483 | |
485 | c *, 5.0350, 5.0330, 4.9640, 4.9410, 4.9050 | |
486 | c *, 4.8714, 4.8408, 4.8482, 4.8518, 4.7611 | |
487 | c *, 4.7157, 4.7049, 4.6834, 4.6939, 4.6206 | |
488 | c *, 4.5931, 4.5723, 4.5498, 4.5210, 4.5337 | |
489 | c *, 4.4726, 4.4898, 4.4251, 4.4013, 4.3784 | |
490 | c *, 4.3514, 4.3264, 4.3530, 4.2846, 4.2567 | |
491 | c *, 4.2339, 4.2171, 4.1919, 4.1757, 4.1503 | |
492 | c *, 4.1329, 4.0924, 4.0797, 4.0624, 4.0446 | |
493 | c *, 4.0210, 4.0010, 3.9812, 3.9627, 3.9402 | |
494 | c *, 3.9221, 3.9029, 3.8925, 3.8680, 3.8439 | |
495 | c *, 3.8244, 3.8094, 3.7816, 3.7719, 3.7396 | |
496 | c *, 3.7271, 3.7011, 3.6849, 3.6804, 3.6640 | |
497 | c *, 3.6496, 3.5954, 3.6042, 3.5697, 3.5552 | |
498 | c *, 3.5463, 3.5265, 3.5212, 3.4992, 3.1432 | |
499 | c *, 3.4630, 3.4483, 3.4484, 3.4240, 3.4127 | |
500 | c *, 3.3595, 3.3767, 3.3582, 3.3059, 3.3258 | |
501 | c *, 3.3182, 3.2582, 3.2740, 3.2621, 3.1984 | |
502 | c *, 3.2408, 3.2043, 3.2087, 3.1652, 3.1720 | |
503 | c *, 3.1476, 3.0754, 3.1236, 3.1494, 3.0653 | |
504 | c *, 3.0757, 3.1030, 3.0692, 3.0466, 3.0392 | |
505 | c *, 3.0346, 2.9933, 3.0146, 2.8774, 2.7972 | |
506 | c *, 2.6463, 2.8068, 2.6869, 2.6643, 2.6673 | |
507 | c *, 2.6684, 2.6409, 2.5520, 2.4341, 2.4978 | |
508 | c *, 2.4534, 2.4417, 2.4293, 2.3283, 2.4978 | |
509 | c *, 2.3042, 2.3111, 2.1996, 2.1996, 2.1960 | |
510 | c *, 2.1469, 2.1223, 2.0791, 2.0599, 1.9891 | |
511 | c *, 2.0027, 1.8678, 1.7850, 1.7868, 1.7362 | |
512 | c *, 1.6478, 1.6254, 1.6166, 1.6478, 1.6516 | |
513 | c *, 1.5554, 1.5727, 1.5656, 1.5290, 1.5149 | |
514 | c *, 1.5554, 1.5615, 1.5522, 1.4863, 1.4745 | |
515 | c *, 1.4745, 1.3991, 1.4139, 1.1686, 1.3458 | |
516 | c *, 1.3587, 1.2425, 1.3171, 1.2828, 1.2153 | |
517 | c *, 1.2679, 1.2766, 1.1410, 1.2514, 1.0852 | |
518 | c *, 1.1445, 1.1672, 1.0624, 1.1089, 1.0899 | |
519 | c *, 1.0171, 1.0399, 1.0645, 1.0417, .9934 | |
520 | c *, 1.0403, 1.0029, 1.0357, 1.0388, 1.0337 | |
521 | c *, 1.0428, 1.0555, 1.0663, 1.0546, 1.0626 | |
522 | c *, 1.0013, 1.0705, 1.0711, 1.0852, 1.0830 | |
523 | c *, 1.1119, 1.0940, 1.0976, 1.0675, 1.1205 | |
524 | c *, 1.0464, 1.0995, .9508, 1.0972, 1.1170 | |
525 | c *, 1.1015, 1.1251, 1.1296, 1.1027, 1.1028 | |
526 | c *, .9271, 1.1363, 1.1120, 1.1455, 1.1192 | |
527 | c *, 1.1498, 1.1491, 1.0108, 1.1282, 1.1550 | |
528 | c *, 1.1617, 1.1354, 1.1586, 1.1631, 1.1378 | |
529 | c *, 1.1656, 1.1684, 1.1389, 1.1694, 1.1690 | |
530 | c *, 1.1702, 1.1705, 1.1390, 1.1714, 1.1716 | |
531 | c *, 1.1326, 1.1517, 1.1697, 1.1731, 1.1716 | |
532 | c *, 1.0867, 1.1537, 1.1698, 1.1642, 1.1638 | |
533 | c *, 1.1199, 1.1635, 1.1713, 1.1631, 1.1601 | |
534 | c *, 1.1340, 1.0823, 1.1598, 1.1754, 1.1572 | |
535 | c *, 1.1565, 1.1567, 1.1631, 1.1538, 1.0852 | |
536 | c *, 1.0342, 1.1645, 1.1452, 1.1099, 1.0940 | |
537 | c *, 1.1185, 1.1470, 1.1388, 1.1424, 1.0705 | |
538 | c *, 1.1374, 1.1199, 1.1262, 1.1142, 1.1205 | |
539 | c *, 1.0867, 1.1095, 1.0734, 1.1234, 1.0925 | |
540 | c *, 1.1135, 1.1076, 1.1070, 1.0955, 1.1027 | |
541 | c *, 1.1073, 1.0630, 1.1007, 1.1185, 1.1133 | |
542 | c *, 1.1128, 1.1028, 1.1033, 1.1086, 1.1120 | |
543 | c *, 1.1013, 1.0991, 1.1086, 1.1028, 1.1027 | |
544 | c *, 1.1064, 1.1036, 1.1135, 1.1139, 1.1136 | |
545 | c *, 1.1145, 1.1166, 1.1152, 1.1168, 1.1172 | |
546 | c *, 1.1230, 1.1187, 1.1254, 1.1224, 1.1329 | |
547 | c *, 1.1216 | |
548 | c */ | |
549 | c data apnecm/ | |
550 | c * 2.1288, 2.2007, 2.2500, 2.3044, 2.3529 | |
551 | c *, 2.9284, 3.5136, 3.6305, 3.7933, 4.1117 | |
552 | c *, 4.5438, 4.9389, 5.1797, 5.3049, 5.4789 | |
553 | c *, 5.6476, 6.2773, 6.9855, 7.6283, 8.2211 | |
554 | c *, 8.7741, 9.2942, 9.7867, 11.5488, 13.7758 | |
555 | c *, 15.0788, 16.8454, 17.9267, 19.4362, 21.2830 | |
556 | c *, 22.9819 | |
557 | c */ | |
558 | c data apnbmx/ | |
559 | c * 1.9462, 1.7481, 1.8881, 1.8019, 1.8627 | |
560 | c *, 1.4745, 1.3231, 1.3762, 1.3351, 1.3505 | |
561 | c *, 1.2841, 1.3086, 1.2565, 1.3038, 1.2218 | |
562 | c *, 1.2952, 1.2127, 1.2035, 1.1968, 1.1951 | |
563 | c *, 1.1604, 1.1794, 1.1747, 1.1724, 1.1595 | |
564 | c *, 1.1637, 1.1477, 1.1520, 1.1510, 1.1478 | |
565 | c *, 1.1466 | |
566 | c */ | |
567 | c data pipecm/ | |
568 | c * 1.1050, 1.1154, 1.1165, 1.1256, 1.1273 | |
569 | c *, 1.1333, 1.1370, 1.1382, 1.1394, 1.1438 | |
570 | c *, 1.1495, 1.1579, 1.1592, 1.1677, 1.1691 | |
571 | c *, 1.1697, 1.1757, 1.1771, 1.1777, 1.1784 | |
572 | c *, 1.1798, 1.1838, 1.1892, 1.1906, 1.1926 | |
573 | c *, 1.1933, 1.1953, 1.1960, 1.1967, 1.1981 | |
574 | c *, 1.1994, 1.2008, 1.2015, 1.2022, 1.2028 | |
575 | c *, 1.2056, 1.2063, 1.2069, 1.2097, 1.2104 | |
576 | c *, 1.2111, 1.2124, 1.2131, 1.2138, 1.2152 | |
577 | c *, 1.2166, 1.2172, 1.2193, 1.2200, 1.2214 | |
578 | c *, 1.2221, 1.2255, 1.2262, 1.2269, 1.2276 | |
579 | c *, 1.2283, 1.2303, 1.2310, 1.2317, 1.2352 | |
580 | c *, 1.2358, 1.2365, 1.2372, 1.2400, 1.2407 | |
581 | c *, 1.2421, 1.2434, 1.2462, 1.2476, 1.2503 | |
582 | c *, 1.2517, 1.2524, 1.2538, 1.2545, 1.2565 | |
583 | c *, 1.2613, 1.2627, 1.2696, 1.2716, 1.2730 | |
584 | c *, 1.2751, 1.2799, 1.2819, 1.2833, 1.2860 | |
585 | c *, 1.2867, 1.2915, 1.2922, 1.2929, 1.2990 | |
586 | c *, 1.3010, 1.3024, 1.3119, 1.3126, 1.3180 | |
587 | c *, 1.3200, 1.3220, 1.3254, 1.3261, 1.3321 | |
588 | c *, 1.3382, 1.3388, 1.3415, 1.3422, 1.3449 | |
589 | c *, 1.3469, 1.3522, 1.3608, 1.3615, 1.3622 | |
590 | c *, 1.3655, 1.3661, 1.3767, 1.3786, 1.3813 | |
591 | c *, 1.3898, 1.3917, 1.3950, 1.4047, 1.4080 | |
592 | c *, 1.4112, 1.4157, 1.4163, 1.4176, 1.4208 | |
593 | c *, 1.4215, 1.4285, 1.4298, 1.4304, 1.4336 | |
594 | c *, 1.4349, 1.4425, 1.4470, 1.4495, 1.4526 | |
595 | c *, 1.4539, 1.4652, 1.4677, 1.4702, 1.4764 | |
596 | c *, 1.4777, 1.4857, 1.4863, 1.4870, 1.4882 | |
597 | c *, 1.4919, 1.4993, 1.4999, 1.5030, 1.5097 | |
598 | c *, 1.5103, 1.5121, 1.5146, 1.5212, 1.5249 | |
599 | c *, 1.5285, 1.5327, 1.5351, 1.5429, 1.5435 | |
600 | c *, 1.5465, 1.5513, 1.5531, 1.5548, 1.5566 | |
601 | c *, 1.5637, 1.5655, 1.5673, 1.5708, 1.5720 | |
602 | c *, 1.5732, 1.5761, 1.5790, 1.5849, 1.5866 | |
603 | c *, 1.5965, 1.5977, 1.5994, 1.6023, 1.6063 | |
604 | c *, 1.6121, 1.6133, 1.6144, 1.6264, 1.6287 | |
605 | c *, 1.6327, 1.6344, 1.6406, 1.6417, 1.6429 | |
606 | c *, 1.6502, 1.6519, 1.6536, 1.6625, 1.6664 | |
607 | c *, 1.6670, 1.6676, 1.6687, 1.6692, 1.6714 | |
608 | c *, 1.6792, 1.6825, 1.6853, 1.6913, 1.6935 | |
609 | c *, 1.6941, 1.6963, 1.6985, 1.6990, 1.7012 | |
610 | c *, 1.7170, 1.7175, 1.7208, 1.7240, 1.7267 | |
611 | c *, 1.7391, 1.7417, 1.7423, 1.7428, 1.7503 | |
612 | c *, 1.7535, 1.7609, 1.7635, 1.7651, 1.7714 | |
613 | c *, 1.7735, 1.7767, 1.7793, 1.7798, 1.7908 | |
614 | c *, 1.7923, 1.7949, 1.7960, 1.8032, 1.8063 | |
615 | c *, 1.8099, 1.8197, 1.8207, 1.8238, 1.8253 | |
616 | c *, 1.8289, 1.8304, 1.8320, 1.8386, 1.8436 | |
617 | c *, 1.8497, 1.8507, 1.8537, 1.8643, 1.8673 | |
618 | c *, 1.8678, 1.8703, 1.8713, 1.8738, 1.8777 | |
619 | c *, 1.8792, 1.8812, 1.8827, 1.8881, 1.8911 | |
620 | c *, 1.8916, 1.8936, 1.8980, 1.9039, 1.9073 | |
621 | c *, 1.9108, 1.9151, 1.9161, 1.9181, 1.9186 | |
622 | c *, 1.9249, 1.9283, 1.9317, 1.9331, 1.9350 | |
623 | c *, 1.9384, 1.9423, 1.9519, 1.9571, 1.9614 | |
624 | c *, 1.9652, 1.9681, 1.9756, 1.9780, 1.9804 | |
625 | c *, 1.9893, 1.9921, 1.9968, 1.9987, 2.0034 | |
626 | c *, 2.0085, 2.0196, 2.0344, 2.0358, 2.0385 | |
627 | c *, 2.0545, 2.0636, 2.0654, 2.0704, 2.0816 | |
628 | c *, 2.0839, 2.0933, 2.1106, 2.1151, 2.1257 | |
629 | c *, 2.1279, 2.1388, 2.1545, 2.1606, 2.1761 | |
630 | c *, 2.1770, 2.1804, 2.1817, 2.1920, 2.2022 | |
631 | c *, 2.2060, 2.2187, 2.2480, 2.2646, 2.2934 | |
632 | c *, 2.3056, 2.3339, 2.3499, 2.3535, 2.3538 | |
633 | c *, 2.3737, 2.3894, 2.4050, 2.4128, 2.4283 | |
634 | c *, 2.4398, 2.4513, 2.4703, 2.4892, 2.5080 | |
635 | c *, 2.5266, 2.5451, 2.5561, 2.5634, 2.5816 | |
636 | c *, 2.5997, 2.6069, 2.6177, 2.6355, 2.6568 | |
637 | c *, 2.6708, 2.6987, 2.7057, 2.7195, 2.7401 | |
638 | c *, 2.7469, 2.7606, 2.7741, 2.7977, 2.8010 | |
639 | c *, 2.8077, 2.8343, 2.8409, 2.8672, 2.8769 | |
640 | c *, 2.8997, 2.9093, 2.9350, 2.9414, 2.9636 | |
641 | c *, 2.9731, 3.0013, 3.0045, 3.0107, 3.0355 | |
642 | c *, 3.0570, 3.0662, 3.0876, 3.0967, 3.1268 | |
643 | c *, 3.1328, 3.1566, 3.1862, 3.2067, 3.2155 | |
644 | c *, 3.2445, 3.2733, 3.3018, 3.3216, 3.3329 | |
645 | c *, 3.3609, 3.3887, 3.4163, 3.4190, 3.4327 | |
646 | c *, 3.4436, 3.4707, 3.4869, 3.4976, 3.5244 | |
647 | c *, 3.5509, 3.6033, 3.6550, 3.7059, 3.7462 | |
648 | c *, 3.9247, 3.9887, 3.9981, 4.4001, 4.4341 | |
649 | c *, 4.8193, 4.8387, 5.2120, 5.2246, 5.3889 | |
650 | c *, 5.5535, 5.5603, 5.7265, 5.8880, 5.8912 | |
651 | c *, 5.9671, 6.1984, 6.2421, 6.5084, 6.6369 | |
652 | c *, 6.9138, 7.5617, 8.1584, 8.7144, 8.9794 | |
653 | c *, 9.2369, 9.7314, 9.9412, 10.2020, 10.6517 | |
654 | c *, 11.4987, 13.7295, 15.0339, 16.6334, 16.8018 | |
655 | c *, 17.8835, 18.1439, 21.2400, 22.9386, 24.1342 | |
656 | c *, 25.2733 | |
657 | c */ | |
658 | c data pipbmx/ | |
659 | c * .4424, .5585, .6180, .7485, .7092 | |
660 | c *, .8058, .7777, .9113, .8974, .9934 | |
661 | c *, 1.0896, 1.2653, 1.3078, 1.5076, 1.5472 | |
662 | c *, 1.5656, 1.7019, 1.7504, 1.7934, 1.7626 | |
663 | c *, 1.8168, 1.9706, 2.0163, 2.1140, 2.1749 | |
664 | c *, 2.0576, 2.1851, 2.1148, 2.1924, 2.1851 | |
665 | c *, 2.3104, 2.3480, 2.2883, 2.1924, 2.3602 | |
666 | c *, 2.4570, 2.3262, 2.2708, 2.4722, 2.3296 | |
667 | c *, 2.5124, 2.3194, 2.4476, 2.4033, 2.5497 | |
668 | c *, 2.5313, 2.5514, 2.5181, 2.5125, 2.5156 | |
669 | c *, 2.5105, 2.4398, 2.4201, 2.5074, 2.4978 | |
670 | c *, 2.4463, 2.4863, 2.4856, 2.4469, 2.5231 | |
671 | c *, 2.3534, 2.5357, 2.4482, 2.3796, 2.3823 | |
672 | c *, 2.3582, 2.3650, 2.3870, 2.1705, 2.1185 | |
673 | c *, 2.2341, 2.2057, 2.2284, 2.1178, 2.1705 | |
674 | c *, 2.0622, 2.0490, 2.0122, 1.9041, 1.8873 | |
675 | c *, 1.9091, 1.7897, 1.7499, 1.7535, 1.8797 | |
676 | c *, 1.8455, 1.6087, 1.6381, 1.6361, 1.6737 | |
677 | c *, 1.5329, 1.5327, 1.4879, 1.4351, 1.4461 | |
678 | c *, 1.3854, 1.3517, 1.3334, 1.3695, 1.2741 | |
679 | c *, 1.2989, 1.2183, 1.2361, 1.2074, 1.1456 | |
680 | c *, 1.1997, 1.1439, 1.1381, 1.0786, 1.1185 | |
681 | c *, 1.0693, 1.1339, .9958, 1.0077, 1.0491 | |
682 | c *, .9997, .9611, .9452, .9141, .9190 | |
683 | c *, .8806, .8885, .9200, .8618, .9680 | |
684 | c *, .8702, .8643, .8263, .8234, .7150 | |
685 | c *, .8412, .8101, .8704, .7767, .8766 | |
686 | c *, .7722, .7694, .7956, .7590, .7128 | |
687 | c *, .7159, .6894, .6894, .7139, .6898 | |
688 | c *, .7412, .6784, .7181, .6947, .7695 | |
689 | c *, .6704, .7110, .6910, .6699, .7199 | |
690 | c *, .6935, .6759, .6843, .6933, .7190 | |
691 | c *, .7159, .6865, .6958, .7017, .6794 | |
692 | c *, .7081, .7004, .7128, .7569, .7269 | |
693 | c *, .7356, .7230, .7332, .7548, .7464 | |
694 | c *, .7695, .7703, .7707, .8288, .8117 | |
695 | c *, .8170, .8054, .7990, .7878, .8556 | |
696 | c *, .8423, .8917, .8664, .8253, .8649 | |
697 | c *, .8627, .8813, .8461, .8938, .9288 | |
698 | c *, .8759, .8345, .8972, .8625, .8818 | |
699 | c *, .8983, .8953, .9080, .9277, .9117 | |
700 | c *, .9080, .8579, .8974, .8970, .9184 | |
701 | c *, .9267, .9184, .9092, .8704, .9146 | |
702 | c *, .9184, .9591, .8649, .9296, .9071 | |
703 | c *, .9449, .9407, .9536, .9721, .9608 | |
704 | c *, .9837, .9233, .9542, 1.0004, .9853 | |
705 | c *, .9934, 1.0005, 1.0061, .9358, 1.0180 | |
706 | c *, 1.0280, 1.0573, .9982, 1.0357, 1.0555 | |
707 | c *, .9877, 1.0682, 1.0694, 1.0817, 1.0720 | |
708 | c *, 1.1046, 1.1041, 1.0596, 1.1013, 1.1217 | |
709 | c *, 1.1217, 1.1291, 1.1113, 1.1270, 1.1354 | |
710 | c *, 1.0799, 1.1199, 1.1439, 1.1459, 1.1213 | |
711 | c *, 1.0907, 1.1368, 1.1488, 1.1008, 1.1500 | |
712 | c *, 1.1156, 1.1502, 1.0947, 1.1410, 1.1427 | |
713 | c *, 1.1480, 1.0729, 1.1634, 1.1319, 1.1480 | |
714 | c *, 1.1295, 1.1213, 1.0783, 1.1044, 1.1027 | |
715 | c *, 1.1217, 1.0998, 1.0418, 1.1031, 1.0712 | |
716 | c *, 1.0720, 1.0600, 1.0540, 1.1180, 1.0490 | |
717 | c *, 1.0464, 1.0380, 1.0124, 1.0202, .9821 | |
718 | c *, 1.0013, .9945, .9997, 1.0187, .9805 | |
719 | c *, .9873, .9796, .9608, 1.0045, .9756 | |
720 | c *, .9695, .9674, .9646, .9662, .9659 | |
721 | c *, .9654, .9491, .9643, .9566, .9654 | |
722 | c *, .9636, .9619, .9676, .9774, .9747 | |
723 | c *, .9916, .9843, .9950, .9932, .9641 | |
724 | c *, .9895, .9861, .9871, .9905, .9916 | |
725 | c *, .9657, .9872, .9770, .9796, .9751 | |
726 | c *, .9637, .9689, .9633, .9635, .9601 | |
727 | c *, .9580, .9707, .9581, .9540, .9505 | |
728 | c *, .9509, .9422, .9494, .9657, .9478 | |
729 | c *, .9513, .9541, .9472, .9439, .9640 | |
730 | c *, .9452, .9397, .9440, .9390, .9416 | |
731 | c *, .9424, .9394, .9334, .9366, .9348 | |
732 | c *, .9338, .9414, .9312, .9115, .9286 | |
733 | c *, .9373, .9287, .9190, .9248, .9232 | |
734 | c *, .9339, .9214, .9201, .9183, .9181 | |
735 | c *, .9170, .9150, .9138, .9071, .9117 | |
736 | c *, .9106, .9092, .9084, .9253, .9053 | |
737 | c *, .9061, .9052, .9132, .9032, .9029 | |
738 | c *, .9008, .8973, .8954, .8928, .9115 | |
739 | c *, .9021, .8938, .8997, .8907, .8921 | |
740 | c *, .8834, .8831, .8795, .8774, .8755 | |
741 | c *, .8745, .8682, .8849, .8649, .8705 | |
742 | c *, .8616, .8684, .8691, .8634, .8687 | |
743 | c *, .8636, .8616, .8581, .8571, .8575 | |
744 | c *, .8582, .8551, .8575, .8582, .8618 | |
745 | c *, .8597, .8619, .8634, .8693, .8645 | |
746 | c *, .8682, .8538, .8759, .8818, .8831 | |
747 | c *, .8853/ | |
748 | c data (pimecm(i),i=1,400)/ | |
749 | c * 1.1046, 1.1133, 1.1394, 1.1425, 1.1495 | |
750 | c *, 1.1579, 1.1585, 1.1592, 1.1598, 1.1677 | |
751 | c *, 1.1691, 1.1731, 1.1777, 1.1784, 1.1798 | |
752 | c *, 1.1831, 1.1858, 1.1879, 1.1892, 1.1906 | |
753 | c *, 1.1940, 1.1967, 1.1994, 1.2008, 1.2015 | |
754 | c *, 1.2028, 1.2069, 1.2076, 1.2090, 1.2097 | |
755 | c *, 1.2111, 1.2124, 1.2131, 1.2159, 1.2166 | |
756 | c *, 1.2179, 1.2186, 1.2200, 1.2214, 1.2234 | |
757 | c *, 1.2269, 1.2283, 1.2296, 1.2303, 1.2317 | |
758 | c *, 1.2324, 1.2352, 1.2358, 1.2365, 1.2372 | |
759 | c *, 1.2407, 1.2421, 1.2441, 1.2462, 1.2476 | |
760 | c *, 1.2510, 1.2517, 1.2524, 1.2545, 1.2551 | |
761 | c *, 1.2579, 1.2586, 1.2593, 1.2607, 1.2613 | |
762 | c *, 1.2627, 1.2634, 1.2641, 1.2662, 1.2668 | |
763 | c *, 1.2682, 1.2696, 1.2710, 1.2716, 1.2730 | |
764 | c *, 1.2758, 1.2778, 1.2806, 1.2819, 1.2826 | |
765 | c *, 1.2833, 1.2847, 1.2888, 1.2915, 1.2922 | |
766 | c *, 1.2929, 1.2963, 1.3004, 1.3024, 1.3038 | |
767 | c *, 1.3058, 1.3065, 1.3072, 1.3112, 1.3119 | |
768 | c *, 1.3126, 1.3146, 1.3180, 1.3187, 1.3200 | |
769 | c *, 1.3207, 1.3220, 1.3227, 1.3254, 1.3261 | |
770 | c *, 1.3301, 1.3308, 1.3321, 1.3335, 1.3362 | |
771 | c *, 1.3375, 1.3388, 1.3415, 1.3422, 1.3455 | |
772 | c *, 1.3482, 1.3515, 1.3522, 1.3555, 1.3562 | |
773 | c *, 1.3615, 1.3628, 1.3641, 1.3655, 1.3681 | |
774 | c *, 1.3694, 1.3760, 1.3773, 1.3786, 1.3898 | |
775 | c *, 1.3917, 1.3963, 1.3995, 1.4002, 1.4015 | |
776 | c *, 1.4047, 1.4099, 1.4112, 1.4163, 1.4176 | |
777 | c *, 1.4183, 1.4196, 1.4208, 1.4228, 1.4234 | |
778 | c *, 1.4272, 1.4292, 1.4298, 1.4304, 1.4393 | |
779 | c *, 1.4425, 1.4470, 1.4476, 1.4489, 1.4507 | |
780 | c *, 1.4539, 1.4589, 1.4596, 1.4608, 1.4652 | |
781 | c *, 1.4658, 1.4683, 1.4727, 1.4739, 1.4758 | |
782 | c *, 1.4764, 1.4777, 1.4795, 1.4808, 1.4833 | |
783 | c *, 1.4851, 1.4857, 1.4876, 1.4882, 1.4894 | |
784 | c *, 1.4901, 1.4919, 1.4925, 1.4938, 1.4968 | |
785 | c *, 1.4974, 1.4981, 1.4993, 1.5030, 1.5054 | |
786 | c *, 1.5060, 1.5066, 1.5072, 1.5097, 1.5103 | |
787 | c *, 1.5109, 1.5146, 1.5152, 1.5182, 1.5206 | |
788 | c *, 1.5212, 1.5224, 1.5231, 1.5249, 1.5309 | |
789 | c *, 1.5327, 1.5345, 1.5357, 1.5387, 1.5405 | |
790 | c *, 1.5411, 1.5429, 1.5435, 1.5441, 1.5465 | |
791 | c *, 1.5489, 1.5507, 1.5513, 1.5519, 1.5548 | |
792 | c *, 1.5590, 1.5614, 1.5637, 1.5655, 1.5673 | |
793 | c *, 1.5702, 1.5708, 1.5761, 1.5802, 1.5825 | |
794 | c *, 1.5831, 1.5843, 1.5849, 1.5866, 1.5890 | |
795 | c *, 1.5936, 1.5948, 1.5977, 1.5994, 1.6017 | |
796 | c *, 1.6052, 1.6058, 1.6087, 1.6092, 1.6138 | |
797 | c *, 1.6173, 1.6178, 1.6190, 1.6236, 1.6258 | |
798 | c *, 1.6276, 1.6287, 1.6298, 1.6315, 1.6327 | |
799 | c *, 1.6367, 1.6372, 1.6378, 1.6400, 1.6406 | |
800 | c *, 1.6423, 1.6480, 1.6513, 1.6519, 1.6525 | |
801 | c *, 1.6547, 1.6553, 1.6614, 1.6631, 1.6637 | |
802 | c *, 1.6648, 1.6653, 1.6659, 1.6681, 1.6692 | |
803 | c *, 1.6703, 1.6720, 1.6726, 1.6731, 1.6742 | |
804 | c *, 1.6792, 1.6798, 1.6803, 1.6820, 1.6825 | |
805 | c *, 1.6853, 1.6858, 1.6864, 1.6875, 1.6935 | |
806 | c *, 1.6957, 1.6963, 1.6979, 1.6990, 1.7067 | |
807 | c *, 1.7072, 1.7121, 1.7127, 1.7132, 1.7186 | |
808 | c *, 1.7202, 1.7208, 1.7213, 1.7224, 1.7235 | |
809 | c *, 1.7240, 1.7267, 1.7278, 1.7294, 1.7342 | |
810 | c *, 1.7348, 1.7391, 1.7407, 1.7428, 1.7476 | |
811 | c *, 1.7503, 1.7535, 1.7545, 1.7609, 1.7614 | |
812 | c *, 1.7667, 1.7714, 1.7735, 1.7741, 1.7767 | |
813 | c *, 1.7772, 1.7798, 1.7824, 1.7835, 1.7871 | |
814 | c *, 1.7929, 1.7960, 1.8001, 1.8032, 1.8058 | |
815 | c *, 1.8130, 1.8182, 1.8212, 1.8218, 1.8238 | |
816 | c *, 1.8253, 1.8258, 1.8289, 1.8320, 1.8370 | |
817 | c *, 1.8386, 1.8436, 1.8446, 1.8507, 1.8512 | |
818 | c *, 1.8557, 1.8588, 1.8638, 1.8678, 1.8713 | |
819 | c *, 1.8762, 1.8777, 1.8792, 1.8807, 1.8827 | |
820 | c *, 1.8857, 1.8881, 1.8886, 1.8936, 1.8975 | |
821 | c *, 1.9010, 1.9063, 1.9132, 1.9186, 1.9205 | |
822 | c *, 1.9249, 1.9254, 1.9273, 1.9283, 1.9317 | |
823 | c *, 1.9326, 1.9374, 1.9423, 1.9451, 1.9495 | |
824 | c *, 1.9519, 1.9538, 1.9614, 1.9709, 1.9733 | |
825 | c *, 1.9747, 1.9799, 1.9884, 1.9903, 1.9987 | |
826 | c *, 2.0034, 2.0117, 2.0164, 2.0182, 2.0196 | |
827 | c *, 2.0330, 2.0335, 2.0358, 2.0540, 2.0581 | |
828 | c *, 2.0636, 2.0654, 2.0798, 2.0812, 2.0816 | |
829 | c */ | |
830 | c data (pimecm(i),i=401,578)/ | |
831 | c * 2.0933, 2.1013, 2.1075, 2.1221, 2.1261 | |
832 | c *, 2.1305, 2.1375, 2.1445, 2.1588, 2.1658 | |
833 | c *, 2.1740, 2.1804, 2.1877, 2.2026, 2.2119 | |
834 | c *, 2.2166, 2.2187, 2.2305, 2.2510, 2.2626 | |
835 | c *, 2.2688, 2.2717, 2.3040, 2.3121, 2.3315 | |
836 | c *, 2.3483, 2.3538, 2.3737, 2.3744, 2.3925 | |
837 | c *, 2.4050, 2.4105, 2.4128, 2.4267, 2.4302 | |
838 | c *, 2.4513, 2.4627, 2.4892, 2.5069, 2.5266 | |
839 | c *, 2.5561, 2.5634, 2.5802, 2.5925, 2.5997 | |
840 | c *, 2.6284, 2.6355, 2.6557, 2.6708, 2.6987 | |
841 | c *, 2.7057, 2.7401, 2.7673, 2.7741, 2.7966 | |
842 | c *, 2.8077, 2.8343, 2.8409, 2.8672, 2.8769 | |
843 | c *, 2.8997, 2.9093, 2.9318, 2.9341, 2.9414 | |
844 | c *, 2.9636, 2.9731, 3.0045, 3.0262, 3.0355 | |
845 | c *, 3.0570, 3.0656, 3.0662, 3.0876, 3.0967 | |
846 | c *, 3.1268, 3.1566, 3.1774, 3.1862, 3.2067 | |
847 | c *, 3.2155, 3.2445, 3.2733, 3.3018, 3.3216 | |
848 | c *, 3.3329, 3.3609, 3.3887, 3.4163, 3.4190 | |
849 | c *, 3.4436, 3.4675, 3.4707, 3.4869, 3.4976 | |
850 | c *, 3.5509, 3.6033, 3.6550, 3.6575, 3.7059 | |
851 | c *, 3.7312, 3.7462, 3.8402, 3.8935, 3.9887 | |
852 | c *, 4.0004, 4.3873, 4.4341, 4.6811, 4.7210 | |
853 | c *, 4.7408, 4.8387, 4.8406, 5.0288, 5.0844 | |
854 | c *, 5.2120, 5.2353, 5.3889, 5.4254, 5.5603 | |
855 | c *, 5.6123, 5.7265, 5.7787, 5.8880, 5.9451 | |
856 | c *, 5.9671, 5.9953, 6.0792, 6.1984, 6.2241 | |
857 | c *, 6.3479, 6.5070, 6.6369, 6.8140, 6.9138 | |
858 | c *, 7.0734, 7.2450, 7.3962, 7.5617, 7.7092 | |
859 | c *, 7.9841, 7.9865, 8.1584, 8.1814, 8.3628 | |
860 | c *, 8.4410, 8.6582, 8.6593, 8.7144, 8.9794 | |
861 | c *, 9.2369, 9.2845, 9.4874, 9.7314, 9.7775 | |
862 | c *, 9.9694, 10.2020, 10.3580, 10.4293, 10.5988 | |
863 | c *, 10.6517, 10.8697, 11.0833, 11.4987, 13.7295 | |
864 | c *, 15.0339, 16.6334, 16.8018, 17.8835, 18.1439 | |
865 | c *, 19.3933, 19.6336, 21.2400, 22.9386, 24.1342 | |
866 | c *, 25.2733, 26.0050, 26.3632 | |
867 | c */ | |
868 | c data (pimbmx(i),i=1,400)/ | |
869 | c * .5314, .5314, .6580, .7092, .7506 | |
870 | c *, .8215, .8579, .8406, .8349, .9441 | |
871 | c *, .9657, 1.0363, 1.0686, 1.1085, 1.1185 | |
872 | c *, 1.1899, 1.1658, 1.2218, 1.2322, 1.2653 | |
873 | c *, 1.3273, 1.2374, 1.3611, 1.3716, 1.3267 | |
874 | c *, 1.3986, 1.4150, 1.3399, 1.4701, 1.4516 | |
875 | c *, 1.4723, 1.4955, 1.4494, 1.4161, 1.4558 | |
876 | c *, 1.4625, 1.5006, 1.5128, 1.5160, 1.4584 | |
877 | c *, 1.4991, 1.4791, 1.4217, 1.4799, 1.4799 | |
878 | c *, 1.4588, 1.4844, 1.4172, 1.4273, 1.4439 | |
879 | c *, 1.4166, 1.4127, 1.3739, 1.3624, 1.3971 | |
880 | c *, 1.3297, 1.3231, 1.3285, 1.2878, 1.2952 | |
881 | c *, 1.2679, 1.2641, 1.2976, 1.2386, 1.2463 | |
882 | c *, 1.3025, 1.2489, 1.2239, 1.1902, 1.2114 | |
883 | c *, 1.1955, 1.2083, 1.1658, 1.1686, 1.1603 | |
884 | c *, 1.1424, 1.1185, 1.1160, 1.1013, 1.1070 | |
885 | c *, 1.1041, 1.0823, 1.0645, 1.0779, 1.0295 | |
886 | c *, 1.0311, .9950, 1.0155, .9950, 1.0029 | |
887 | c *, .9767, 1.0249, .9853, .9657, .9906 | |
888 | c *, .9608, .9591, .9805, .9458, .9575 | |
889 | c *, .9422, .9407, .9558, .9473, .9271 | |
890 | c *, .9132, .9236, .9219, .9219, .9575 | |
891 | c *, .9097, .9224, .9094, .9202, .8903 | |
892 | c *, .9224, .9160, .9122, .9094, .8956 | |
893 | c *, .9167, .9271, .9398, .9184, .9277 | |
894 | c *, .9310, .9184, .9354, .9219, .9300 | |
895 | c *, .9303, .9508, .9451, .9680, .9399 | |
896 | c *, .9420, .9387, .9667, .9248, .9915 | |
897 | c *, .9698, .9626, 1.0187, .9659, .9441 | |
898 | c *, .9766, .9694, 1.0061, .9789, .9757 | |
899 | c *, .9803, 1.0495, 1.0218, 1.0026, 1.0706 | |
900 | c *, 1.0277, 1.0302, 1.0690, 1.0552, 1.0406 | |
901 | c *, 1.0801, 1.0522, 1.1284, 1.1226, 1.1270 | |
902 | c *, 1.0934, 1.0968, 1.1381, 1.0915, 1.1543 | |
903 | c *, 1.1775, 1.2101, 1.1247, 1.1506, 1.1781 | |
904 | c *, 1.1546, 1.1968, 1.1604, 1.1697, 1.2072 | |
905 | c *, 1.2058, 1.2489, 1.1912, 1.2141, 1.1928 | |
906 | c *, 1.1986, 1.1978, 1.2256, 1.2565, 1.2029 | |
907 | c *, 1.2127, 1.1997, 1.2035, 1.1791, 1.1853 | |
908 | c *, 1.1914, 1.1927, 1.2616, 1.2303, 1.1848 | |
909 | c *, 1.1507, 1.2399, 1.1510, 1.1898, 1.1499 | |
910 | c *, 1.1337, 1.1433, 1.1089, 1.1517, 1.2008 | |
911 | c *, 1.1598, 1.1030, 1.0798, 1.1165, 1.1213 | |
912 | c *, 1.1190, 1.0861, 1.0776, 1.0653, 1.0730 | |
913 | c *, 1.1106, 1.0908, 1.0945, 1.0650, 1.0882 | |
914 | c *, 1.0720, 1.0540, 1.0925, 1.0718, 1.1377 | |
915 | c *, 1.0867, 1.0802, 1.1090, 1.0936, 1.0915 | |
916 | c *, 1.1316, 1.1575, 1.1202, 1.1275, 1.0908 | |
917 | c *, 1.1142, 1.1507, 1.1890, 1.1863, 1.2114 | |
918 | c *, 1.1142, 1.2001, 1.1312, 1.1930, 1.2563 | |
919 | c *, 1.2489, 1.2919, 1.2652, 1.2425, 1.2149 | |
920 | c *, 1.2386, 1.2616, 1.3207, 1.3299, 1.2816 | |
921 | c *, 1.3296, 1.3562, 1.2957, 1.3831, 1.3028 | |
922 | c *, 1.3806, 1.3704, 1.3202, 1.4082, 1.3482 | |
923 | c *, 1.3730, 1.3351, 1.3886, 1.3348, 1.3957 | |
924 | c *, 1.3929, 1.3790, 1.3791, 1.3946, 1.3641 | |
925 | c *, 1.3877, 1.3566, 1.3564, 1.3762, 1.3253 | |
926 | c *, 1.3739, 1.3252, 1.3234, 1.3423, 1.3159 | |
927 | c *, 1.2643, 1.2700, 1.2584, 1.2466, 1.2841 | |
928 | c *, 1.2191, 1.2052, 1.2794, 1.2835, 1.2374 | |
929 | c *, 1.2061, 1.1936, 1.1947, 1.2476, 1.1481 | |
930 | c *, 1.1496, 1.2187, 1.1427, 1.2101, 1.1192 | |
931 | c *, 1.1445, 1.1371, 1.1613, 1.0994, 1.1020 | |
932 | c *, 1.1233, 1.1041, 1.0896, 1.0890, 1.0880 | |
933 | c *, 1.0802, 1.0946, 1.0940, 1.0670, 1.0761 | |
934 | c *, 1.0817, 1.0702, 1.0723, 1.0779, 1.0789 | |
935 | c *, 1.0733, 1.0771, 1.0767, 1.0710, 1.0705 | |
936 | c *, 1.0720, 1.0729, 1.0472, 1.0959, 1.0788 | |
937 | c *, 1.0739, 1.0645, 1.0681, 1.0749, 1.0789 | |
938 | c *, 1.0799, 1.0795, 1.0763, 1.0588, 1.0805 | |
939 | c *, 1.0724, 1.0660, 1.0663, 1.0801, 1.0988 | |
940 | c *, 1.0797, 1.0816, 1.0687, 1.0690, 1.0724 | |
941 | c *, 1.0721, 1.0777, 1.0650, 1.0715, 1.0617 | |
942 | c *, 1.0690, 1.0678, 1.0714, 1.0323, 1.0672 | |
943 | c *, 1.0630, 1.0564, 1.0600, 1.0636, 1.0499 | |
944 | c *, 1.0479, 1.0481, 1.0463, 1.0510, 1.0450 | |
945 | c *, 1.0171, 1.0517, 1.0498, 1.0510, 1.0412 | |
946 | c *, 1.0374, 1.0501, 1.0388, 1.0449, 1.0449 | |
947 | c *, 1.0311, 1.0510, 1.0505, 1.0449, 1.0567 | |
948 | c *, 1.0572, 1.0495, 1.0632, 1.0461, 1.0403 | |
949 | c */ | |
950 | c data (pimbmx(i),i=401,578)/ | |
951 | c * 1.0629, 1.0642, 1.0525, 1.0709, 1.0585 | |
952 | c *, 1.0612, 1.0717, 1.0731, 1.0660, 1.0761 | |
953 | c *, 1.0696, 1.0767, 1.0755, 1.0720, 1.0615 | |
954 | c *, 1.0665, 1.0660, 1.0714, 1.0670, 1.0650 | |
955 | c *, 1.0627, 1.0621, 1.0499, 1.0499, 1.0457 | |
956 | c *, 1.0405, 1.0372, 1.0412, 1.0331, 1.0299 | |
957 | c *, 1.0299, 1.0241, 1.0321, 1.0306, 1.0232 | |
958 | c *, 1.0252, 1.0171, 1.0211, 1.0226, 1.0178 | |
959 | c *, 1.0083, 1.0157, 1.0138, .9918, 1.0140 | |
960 | c *, 1.0028, 1.0120, 1.0077, 1.0088, .9979 | |
961 | c *, 1.0052, 1.0025, .9876, .9987, .9990 | |
962 | c *, .9944, .9845, .9918, .9777, .9892 | |
963 | c *, .9754, .9856, .9812, .9816, .9831 | |
964 | c *, .9674, .9800, .9781, .9679, .9756 | |
965 | c *, .9805, .9759, .9739, .9672, .9718 | |
966 | c *, .9701, .9682, .9632, .9664, .9538 | |
967 | c *, .9647, .9628, .9613, .9596, .9449 | |
968 | c *, .9585, .9572, .9553, .9540, .9624 | |
969 | c *, .9527, .9468, .9513, .9525, .9503 | |
970 | c *, .9479, .9459, .9438, .9424, .9422 | |
971 | c *, .9407, .9457, .9399, .9385, .9368 | |
972 | c *, .9409, .9248, .9209, .9034, .8974 | |
973 | c *, .9150, .9089, .9145, .9146, .9080 | |
974 | c *, .9044, .9082, .9069, .9062, .8938 | |
975 | c *, .9034, .9045, .9011, .8921, .8979 | |
976 | c *, .8925, .8982, .8975, .8976, .8947 | |
977 | c *, .8962, .8932, .8916, .8913, .8888 | |
978 | c *, .8889, .8892, .8880, .8881, .8842 | |
979 | c *, .8815, .8826, .8832, .8813, .8826 | |
980 | c *, .8791, .8835, .8831, .8833, .8782 | |
981 | c *, .8789, .8813, .8797, .8780, .8744 | |
982 | c *, .8782, .8842, .8777, .8777, .8789 | |
983 | c *, .8815, .8800, .8839, .8740, .8732 | |
984 | c *, .8751, .8784, .8757, .8779, .8630 | |
985 | c *, .8802, .8775, .8851, .8874, .8903 | |
986 | c *, .8935, .8965, .8965 | |
987 | c */ | |
988 | c data kmpecm/ | |
989 | c * 1.4691, 1.4720, 1.4750, 1.4780, 1.4811 | |
990 | c *, 1.4837, 1.4843, 1.4860, 1.4876, 1.4910 | |
991 | c *, 1.4944, 1.4979, 1.5014, 1.5032, 1.5050 | |
992 | c *, 1.5087, 1.5091, 1.5124, 1.5132, 1.5162 | |
993 | c *, 1.5170, 1.5200, 1.5220, 1.5239, 1.5278 | |
994 | c *, 1.5318, 1.5354, 1.5358, 1.5362, 1.5378 | |
995 | c *, 1.5399, 1.5440, 1.5523, 1.5607, 1.5654 | |
996 | c *, 1.5688, 1.5775, 1.5784, 1.5863, 1.5916 | |
997 | c *, 1.5947, 1.6023, 1.6050, 1.6055, 1.6086 | |
998 | c *, 1.6145, 1.6159, 1.6172, 1.6191, 1.6236 | |
999 | c *, 1.6319, 1.6328, 1.6332, 1.6420, 1.6461 | |
1000 | c *, 1.6466, 1.6522, 1.6563, 1.6582, 1.6614 | |
1001 | c *, 1.6642, 1.6694, 1.6712, 1.6717, 1.6768 | |
1002 | c *, 1.6806, 1.6811, 1.6839, 1.6843, 1.6867 | |
1003 | c *, 1.6885, 1.6961, 1.6965, 1.7003, 1.7022 | |
1004 | c *, 1.7083, 1.7087, 1.7172, 1.7177, 1.7181 | |
1005 | c *, 1.7229, 1.7243, 1.7276, 1.7342, 1.7374 | |
1006 | c *, 1.7436, 1.7459, 1.7483, 1.7539, 1.7610 | |
1007 | c *, 1.7629, 1.7633, 1.7718, 1.7770, 1.7789 | |
1008 | c *, 1.7793, 1.7817, 1.7840, 1.7873, 1.7892 | |
1009 | c *, 1.8028, 1.8037, 1.8075, 1.8135, 1.8140 | |
1010 | c *, 1.8219, 1.8247, 1.8261, 1.8308, 1.8369 | |
1011 | c *, 1.8401, 1.8406, 1.8410, 1.8475, 1.8480 | |
1012 | c *, 1.8489, 1.8540, 1.8559, 1.8605, 1.8647 | |
1013 | c *, 1.8721, 1.8744, 1.8767, 1.8785, 1.8836 | |
1014 | c *, 1.8951, 1.8955, 1.8983, 1.9001, 1.9029 | |
1015 | c *, 1.9065, 1.9184, 1.9202, 1.9243, 1.9348 | |
1016 | c *, 1.9393, 1.9411, 1.9434, 1.9483, 1.9547 | |
1017 | c *, 1.9628, 1.9637, 1.9659, 1.9699, 1.9798 | |
1018 | c *, 1.9838, 1.9923, 1.9958, 2.0047, 2.0127 | |
1019 | c *, 2.0149, 2.0162, 2.0255, 2.0277, 2.0308 | |
1020 | c *, 2.0417, 2.0430, 2.0487, 2.0579, 2.0653 | |
1021 | c *, 2.0679, 2.0713, 2.0813, 2.0882, 2.0925 | |
1022 | c *, 2.1028, 2.1105, 2.1126, 2.1211, 2.1232 | |
1023 | c *, 2.1258, 2.1351, 2.1444, 2.1507, 2.1528 | |
1024 | c *, 2.1654, 2.1675, 2.1687, 2.1737, 2.1837 | |
1025 | c *, 2.1862, 2.1937, 2.2044, 2.2131, 2.2143 | |
1026 | c *, 2.2274, 2.2356, 2.2478, 2.2546, 2.2659 | |
1027 | c *, 2.2756, 2.2836, 2.2976, 2.2996, 2.3162 | |
1028 | c *, 2.3166, 2.3296, 2.3335, 2.3363, 2.3535 | |
1029 | c *, 2.3562, 2.3725, 2.3729, 2.3748, 2.3887 | |
1030 | c *, 2.3918, 2.3933, 2.4006, 2.4109, 2.4174 | |
1031 | c *, 2.4223, 2.4299, 2.4352, 2.4488, 2.4518 | |
1032 | c *, 2.4675, 2.4705, 2.4861, 2.4887, 2.4936 | |
1033 | c *, 2.5046, 2.5230, 2.5412, 2.5593, 2.5701 | |
1034 | c *, 2.5773, 2.5952, 2.6023, 2.6059, 2.6130 | |
1035 | c *, 2.6306, 2.6447, 2.6482, 2.6656, 2.6795 | |
1036 | c *, 2.6829, 2.7002, 2.7173, 2.7848, 2.8540 | |
1037 | c *, 2.9216, 2.9407, 2.9878, 3.0096, 3.0250 | |
1038 | c *, 3.0526, 3.1783, 3.2191, 3.2993, 3.3887 | |
1039 | c *, 3.5239, 3.6924, 3.7800, 3.9967, 4.0200 | |
1040 | c *, 4.1348, 4.4617, 4.4826, 4.7663, 4.8636 | |
1041 | c *, 4.9779, 5.1080, 5.1263, 5.2349, 5.3237 | |
1042 | c *, 5.4110, 5.5479, 5.5816, 5.8281, 5.9080 | |
1043 | c *, 6.0801, 6.2173, 6.3664, 6.6545, 6.9306 | |
1044 | c *, 7.2610, 7.5770, 7.7605, 7.8207, 7.9985 | |
1045 | c *, 8.1725, 8.2297, 8.4546, 8.7275, 8.9922 | |
1046 | c *, 9.2493, 9.4994, 9.7431, 9.9809, 10.2131 | |
1047 | c *, 11.5086, 13.7378, 15.0415, 16.8085, 17.8898 | |
1048 | c *, 19.3991, 21.2453, 22.9435, 24.1389 | |
1049 | c */ | |
1050 | c data kmpbmx/ | |
1051 | c * 1.9033, 1.7662, 1.7298, 1.7544, 1.5461 | |
1052 | c *, 1.6991, 1.6205, 1.5898, 1.5817, 1.5023 | |
1053 | c *, 1.5554, 1.5086, 1.5065, 1.4948, 1.4799 | |
1054 | c *, 1.4927, 1.4854, 1.6136, 1.6017, 1.7312 | |
1055 | c *, 1.4884, 1.7075, 1.5574, 1.5260, 1.5002 | |
1056 | c *, 1.4571, 1.3991, 1.3219, 1.3434, 1.3635 | |
1057 | c *, 1.3315, 1.2966, 1.2213, 1.1604, 1.1930 | |
1058 | c *, 1.1382, 1.1354, 1.1583, 1.1185, 1.1156 | |
1059 | c *, 1.0645, 1.0247, 1.0969, 1.0779, 1.1083 | |
1060 | c *, 1.0749, 1.0155, 1.0464, 1.0754, 1.0202 | |
1061 | c *, 1.0470, 1.0615, 1.0540, 1.0357, 1.0605 | |
1062 | c *, 1.0428, 1.0265, 1.0187, 1.0611, 1.0615 | |
1063 | c *, 1.0295, 1.0651, 1.0823, 1.0077, 1.0764 | |
1064 | c *, 1.1298, 1.1013, 1.1080, 1.1070, 1.1368 | |
1065 | c *, 1.1382, 1.1410, 1.1418, 1.1312, 1.1466 | |
1066 | c *, 1.1378, 1.1241, 1.0896, 1.1264, 1.0720 | |
1067 | c *, 1.1368, 1.1013, 1.1442, 1.1466, 1.1612 | |
1068 | c *, 1.1726, 1.1755, 1.1576, 1.1726, 1.1848 | |
1069 | c *, 1.1617, 1.2029, 1.1686, 1.2274, 1.2252 | |
1070 | c *, 1.1726, 1.2256, 1.2244, 1.2008, 1.2332 | |
1071 | c *, 1.2828, 1.2548, 1.2412, 1.2887, 1.2853 | |
1072 | c *, 1.2653, 1.2540, 1.2527, 1.2449, 1.2118 | |
1073 | c *, 1.1781, 1.1902, 1.1767, 1.1859, 1.1410 | |
1074 | c *, 1.1594, 1.1800, 1.1185, 1.1256, 1.1354 | |
1075 | c *, 1.1131, 1.1562, 1.1143, 1.1382, 1.0841 | |
1076 | c *, 1.0587, 1.0632, 1.0311, 1.0617, 1.0150 | |
1077 | c *, 1.0309, 1.0174, 1.0110, 1.0171, .9961 | |
1078 | c *, .9719, .9938, .9869, .9953, .9966 | |
1079 | c *, 1.0189, .9977, .9918, .9948, 1.0034 | |
1080 | c *, .9737, 1.0066, 1.0137, 1.0041, 1.0171 | |
1081 | c *, 1.0303, 1.0223, 1.0322, .9751, 1.0331 | |
1082 | c *, 1.0118, 1.0403, 1.0399, 1.0429, 1.0171 | |
1083 | c *, 1.0279, 1.0467, 1.0414, 1.0217, 1.0432 | |
1084 | c *, 1.0379, 1.0280, 1.0351, 1.0171, 1.0278 | |
1085 | c *, 1.0314, 1.0213, 1.0133, 1.0240, 1.0080 | |
1086 | c *, 1.0018, .9964, .9984, .9889, .9910 | |
1087 | c *, .9903, .9801, .9852, .9853, .9847 | |
1088 | c *, .9800, .9832, .9770, .9749, .9754 | |
1089 | c *, .9723, .9741, .9744, .9738, .9751 | |
1090 | c *, .9780, .9738, .9738, .9684, .9712 | |
1091 | c *, .9669, .9680, .9671, .9576, .9619 | |
1092 | c *, .9624, .9585, .9588, .9586, .9525 | |
1093 | c *, .9253, .9518, .9503, .9491, .9463 | |
1094 | c *, .9476, .9420, .9458, .9434, .9341 | |
1095 | c *, .9444, .9412, .9393, .9395, .9210 | |
1096 | c *, .9370, .9358, .8974, .9400, .9342 | |
1097 | c *, .9305, .9141, .9271, .9267, .9228 | |
1098 | c *, .9233, .9219, .9247, .9296, .9253 | |
1099 | c *, .9089, .8992, .8946, .8992, .9062 | |
1100 | c *, .9069, .8874, .9219, .8746, .8795 | |
1101 | c *, .8740, .8704, .8785, .8795, .8667 | |
1102 | c *, .8849, .8510, .8463, .8612, .8292 | |
1103 | c *, .8292, .8387, .8273, .8273, .8292 | |
1104 | c *, .8292, .8349, .8234, .8349, .8176 | |
1105 | c *, .8292, .8180, .8193, .8154, .8130 | |
1106 | c *, .8121, .8145, .8078, .7959, .8088 | |
1107 | c *, .8075, .8064, .8056, .8086, .8048 | |
1108 | c *, .8080, .8068, .8050, .8042, .8050 | |
1109 | c *, .8054, .8064, .8096, .8095, .8107 | |
1110 | c *, .8135, .8234, .8238, .8263 | |
1111 | c */ | |
1112 | c data kmnecm/ | |
1113 | c * 1.6212, 1.6781, 1.7053, 1.7863, 1.7882 | |
1114 | c *, 1.8271, 1.9025, 2.0152, 2.1237, 2.2157 | |
1115 | c *, 2.4252, 2.6054, 2.6160, 2.9441, 3.5279 | |
1116 | c *, 3.6965, 4.0245, 4.4666, 4.8690, 5.1136 | |
1117 | c *, 5.2406, 5.4169, 5.5877, 5.9144, 6.2241 | |
1118 | c *, 6.9381, 7.5852, 8.1813, 8.7369, 9.2592 | |
1119 | c *, 9.7536, 10.2241, 11.5209, 13.7525, 15.0575 | |
1120 | c *, 16.8264, 17.9089, 19.4198, 21.2680, 22.9680 | |
1121 | c *, 24.1646 | |
1122 | c */ | |
1123 | c data kmnbmx/ | |
1124 | c * .9132, .9966, .9772, .9966, 1.1382 | |
1125 | c *, 1.0794, .9674, .9167, .8795, .8500 | |
1126 | c *, .8482, .8444, .8330, .8078, .8349 | |
1127 | c *, .8368, .7919, .8146, .8019, .8156 | |
1128 | c *, .7999, .7961, .8038, .8038, .7949 | |
1129 | c *, .7868, .7864, .7910, .7850, .7870 | |
1130 | c *, .7906, .7885, .7945, .7971, .8011 | |
1131 | c *, .8003, .8029, .8082, .8088, .8156 | |
1132 | c *, .8189 | |
1133 | c */ | |
1134 | c data kppecm/ | |
1135 | c * 1.4447, 1.4516, 1.4522, 1.4585, 1.4663 | |
1136 | c *, 1.4750, 1.5036, 1.5050, 1.5091, 1.5239 | |
1137 | c *, 1.5378, 1.5523, 1.5654, 1.5784, 1.5916 | |
1138 | c *, 1.5929, 1.6014, 1.6037, 1.6050, 1.6150 | |
1139 | c *, 1.6159, 1.6191, 1.6259, 1.6264, 1.6268 | |
1140 | c *, 1.6328, 1.6378, 1.6461, 1.6517, 1.6587 | |
1141 | c *, 1.6605, 1.6652, 1.6792, 1.6843, 1.6853 | |
1142 | c *, 1.6928, 1.7073, 1.7101, 1.7210, 1.7294 | |
1143 | c *, 1.7374, 1.7422, 1.7483, 1.7539, 1.7643 | |
1144 | c *, 1.7662, 1.7704, 1.7789, 1.7793, 1.7864 | |
1145 | c *, 1.7897, 1.8028, 1.8070, 1.8135, 1.8191 | |
1146 | c *, 1.8327, 1.8355, 1.8373, 1.8517, 1.8587 | |
1147 | c *, 1.8605, 1.8679, 1.8725, 1.8813, 1.8836 | |
1148 | c *, 1.9038, 1.9070, 1.9289, 1.9298, 1.9321 | |
1149 | c *, 1.9524, 1.9533, 1.9749, 1.9807, 1.9950 | |
1150 | c *, 1.9972, 1.9994, 2.0074, 2.0193, 2.0435 | |
1151 | c *, 2.0492, 2.0635, 2.0653, 2.0851, 2.1040 | |
1152 | c *, 2.1066, 2.1083, 2.1279, 2.1296, 2.1490 | |
1153 | c *, 2.1507, 2.1717, 2.1908, 2.1924, 2.2110 | |
1154 | c *, 2.2131, 2.2213, 2.2319, 2.2335, 2.2538 | |
1155 | c *, 2.2724, 2.2740, 2.2940, 2.3123, 2.3138 | |
1156 | c *, 2.3375, 2.3531, 2.3725, 2.3903, 2.3918 | |
1157 | c *, 2.4109, 2.4197, 2.4299, 2.4413, 2.4488 | |
1158 | c *, 2.4675, 2.4861, 2.5046, 2.5230, 2.5266 | |
1159 | c *, 2.5412, 2.5521, 2.5593, 2.5773, 2.5952 | |
1160 | c *, 2.6130, 2.6306, 2.6482, 2.6656, 2.6829 | |
1161 | c *, 2.7002, 2.7173, 3.0096, 3.0556, 3.1754 | |
1162 | c *, 3.3887, 3.5239, 3.7800, 4.0200, 4.1348 | |
1163 | c *, 4.4617, 4.6469, 4.7663, 4.8636, 4.9591 | |
1164 | c *, 5.1263, 5.2349, 5.4110, 5.5816, 5.7308 | |
1165 | c *, 5.9080, 6.0646, 6.2173, 6.9306, 7.5770 | |
1166 | c *, 7.8207, 8.1725, 8.7275, 9.2493, 9.7431 | |
1167 | c *, 10.2131, 11.5086, 13.7378, 15.0415, 16.6402 | |
1168 | c *, 16.8085, 17.8898, 18.1501, 19.3991, 21.2453 | |
1169 | c *, 22.9435, 24.1389 | |
1170 | c */ | |
1171 | c data kppbmx/ | |
1172 | c * .5412, .6308, .6024, .6050, .5971 | |
1173 | c *, .6037, .6232, .6103, .6482, .6601 | |
1174 | c *, .6386, .6466, .6438, .6204, .6482 | |
1175 | c *, .6358, .6333, .6445, .6443, .6346 | |
1176 | c *, .6410, .6227, .6283, .6308, .6403 | |
1177 | c *, .6290, .6457, .5984, .6333, .5955 | |
1178 | c *, .5944, .6295, .6346, .6090, .6433 | |
1179 | c *, .6383, .6482, .6425, .6543, .6588 | |
1180 | c *, .6652, .6768, .6730, .6723, .6815 | |
1181 | c *, .7040, .6898, .7014, .7001, .7181 | |
1182 | c *, .7130, .7159, .7067, .7440, .7345 | |
1183 | c *, .7365, .7485, .7382, .7474, .7574 | |
1184 | c *, .7588, .7559, .7590, .7582, .7668 | |
1185 | c *, .7592, .7682, .7661, .7697, .7548 | |
1186 | c *, .7661, .7626, .7626, .7563, .7590 | |
1187 | c *, .7578, .7611, .7557, .7555, .7506 | |
1188 | c *, .7498, .7517, .7508, .7540, .7464 | |
1189 | c *, .7538, .7512, .7527, .7534, .7527 | |
1190 | c *, .7565, .7521, .7529, .7525, .7444 | |
1191 | c *, .7517, .7334, .7485, .7491, .7510 | |
1192 | c *, .7466, .7476, .7478, .7472, .7485 | |
1193 | c *, .7378, .7451, .7468, .7474, .7476 | |
1194 | c *, .7459, .7410, .7461, .7457, .7414 | |
1195 | c *, .7464, .7457, .7444, .7444, .7444 | |
1196 | c *, .7442, .7291, .7421, .7429, .7421 | |
1197 | c *, .7397, .7386, .7373, .7389, .7384 | |
1198 | c *, .7384, .7386, .7378, .7569, .7878 | |
1199 | c *, .7548, .7356, .7526, .7421, .7715 | |
1200 | c *, .7568, .7590, .7777, .7421, .7632 | |
1201 | c *, .7464, .7442, .7548, .7367, .7736 | |
1202 | c *, .7378, .7421, .7455, .7502, .7510 | |
1203 | c *, .7653, .7529, .7580, .7544, .7614 | |
1204 | c *, .7605, .7661, .7703, .7805, .7883 | |
1205 | c *, .7850, .7907, .7611, .7961, .8023 | |
1206 | c *, .8068, .8111 | |
1207 | c */ | |
1208 | c data kpnecm/ | |
1209 | c * 1.6875, 1.7430, 1.7670, 1.7816, 1.7905 | |
1210 | c *, 1.8144, 1.8383, 1.8615, 1.8749, 1.8846 | |
1211 | c *, 1.9080, 1.9308, 1.9345, 1.9535, 1.9760 | |
1212 | c *, 1.9974, 1.9983, 2.0205, 2.0460, 2.0647 | |
1213 | c *, 2.0678, 2.0864, 2.1066, 2.1079, 2.1109 | |
1214 | c *, 2.1292, 2.1322, 2.1504, 2.1533, 2.1743 | |
1215 | c *, 2.1922, 2.1951, 2.2157, 2.2239, 2.2334 | |
1216 | c *, 2.2362, 2.2565, 2.2739, 2.2767, 2.2967 | |
1217 | c *, 2.3138, 2.3166, 2.3402, 2.3559, 2.3753 | |
1218 | c *, 2.3919, 2.3946, 2.4138, 2.4328, 2.4517 | |
1219 | c *, 2.4704, 2.4891, 2.5076, 2.5259, 2.5442 | |
1220 | c *, 2.5551, 2.5623, 2.5803, 2.5982, 2.6160 | |
1221 | c *, 2.6337, 2.6513, 2.6687, 2.6861, 2.7033 | |
1222 | c *, 2.7204, 3.5279, 4.0245, 4.4666, 4.8690 | |
1223 | c *, 5.2406, 5.4169, 5.5877, 5.9144, 6.2241 | |
1224 | c *, 6.9381, 7.5852, 8.1813, 8.7369, 9.2592 | |
1225 | c *, 9.7536, 10.2241, 11.5209, 13.7525, 15.0575 | |
1226 | c *, 16.8264, 17.9089, 19.4198, 21.2680, 22.9680 | |
1227 | c *, 24.1646 | |
1228 | c */ | |
1229 | c data kpnbmx/ | |
1230 | c * .7024, .7324, .7485, .7527, .7680 | |
1231 | c *, .7758, .8100, .8224, .7611, .8151 | |
1232 | c *, .8031, .7915, .7674, .7842, .7822 | |
1233 | c *, .7590, .7791, .7767, .7758, .7734 | |
1234 | c *, .7754, .7709, .7674, .7713, .7742 | |
1235 | c *, .7752, .7748, .7721, .7680, .7707 | |
1236 | c *, .7674, .7713, .7715, .7695, .7684 | |
1237 | c *, .7734, .7682, .7709, .7672, .7659 | |
1238 | c *, .7653, .7653, .7506, .7626, .7624 | |
1239 | c *, .7701, .7588, .7622, .7592, .7491 | |
1240 | c *, .7588, .7574, .7592, .7582, .7571 | |
1241 | c *, .7464, .7559, .7538, .7529, .7529 | |
1242 | c *, .7534, .7538, .7487, .7487, .7498 | |
1243 | c *, .7474, .7464, .7485, .7464, .7485 | |
1244 | c *, .7464, .7565, .7442, .7485, .7545 | |
1245 | c *, .7555, .7542, .7634, .7653, .7686 | |
1246 | c *, .7671, .7723, .7695, .7785, .7824 | |
1247 | c *, .7905, .7927, .7923, .8052, .8100 | |
1248 | c *, .8137 | |
1249 | c */ | |
1250 | cc data lpecm/ | |
1251 | cc * 2.0577, 2.0583, 2.0593, 2.0595, 2.0609 | |
1252 | cc *, 2.0617, 2.0629, 2.0642, 2.0643, 2.0647 | |
1253 | cc *, 2.0666, 2.0671, 2.0683, 2.0692, 2.0709 | |
1254 | cc *, 2.0720, 2.0783, 2.0818, 2.0935, 2.1022 | |
1255 | cc *, 2.1058, 2.1117, 2.1326, 2.1559, 2.1731 | |
1256 | cc *, 2.1811, 2.2079, 2.2360, 2.2505, 2.3107 | |
1257 | cc *, 2.3733, 2.4534, 2.4695, 2.9278, 5.2476 | |
1258 | cc */ | |
1259 | cc data lpbmx/ | |
1260 | cc * 2.5793, 2.3937, 1.9381, 2.3736, 2.0342 | |
1261 | cc *, 2.2068, 1.9381, 1.8797, 1.7841, 1.7930 | |
1262 | cc *, 1.2676, 1.6641, 1.4604, 1.0403, 1.3470 | |
1263 | cc *, 1.0420, .8722, .8368, .5323, .5352 | |
1264 | cc *, .8867, .5382, .8106, .8043, .8058 | |
1265 | cc *, .8630, .7031, .7777, .9441, 1.0403 | |
1266 | cc *, .9772, .6628, 1.1968, 1.0555, 1.0495 | |
1267 | cc */ | |
1268 | c data pi1ecm/ | |
1269 | c * .45000, .55000, .62500, .68750, .73750 | |
1270 | c *, .78750, .83750, .88750, .95000, 1.02500 | |
1271 | c *, 1.10000, 1.17500 | |
1272 | c */ | |
1273 | c data pi1bmx/ | |
1274 | c * .48533, .47873, .55852, .58087, .55279 | |
1275 | c *, .49185, .47203, .46181, .46181, .47203 | |
1276 | c *, .41459, .40684 | |
1277 | c */ | |
1278 | c data pi2ecm/ | |
1279 | c * .45000, .55000, .62500, .68750, .73750 | |
1280 | c *, .78750, .83750, .88750, .95000, 1.02500 | |
1281 | c *, 1.10000, 1.17500 | |
1282 | c */ | |
1283 | c data pi2bmx/ | |
1284 | c * .45135, .80385, .87767, 1.28779, 1.81771 | |
1285 | c *, 2.07296, 1.55536, .95912, .94407, .77768 | |
1286 | c *, .70467, .54408 | |
1287 | c */ | |
1288 | c data pi3ecm/ | |
1289 | c * .45000, .55000, .61250, .65000, .67500 | |
1290 | c *, .70000, .72500, .75000, .77500, .80000 | |
1291 | c *, .82500, .85000, .87500, .90000, .93750 | |
1292 | c *, 1.02500, 1.10000, 1.17500 | |
1293 | c */ | |
1294 | c data pi3bmx/ | |
1295 | c * .50463, .79988, 1.13541, 1.30497, 1.34462 | |
1296 | c *, 1.52957, 1.63809, 1.94216, 1.98511, 1.66316 | |
1297 | c *, 1.69163, 1.62933, 1.06152, .96574, .83302 | |
1298 | c *, .71365, .66517, .59173 | |
1299 | c */ | |
1300 | c data pi4ecm/ | |
1301 | c * .28500, .29500, .30500, .31500, .37000 | |
1302 | c *, .46000, .54375, .61250, .65000, .67500 | |
1303 | c *, .70000, .72500, .75000, .77500, .80000 | |
1304 | c *, .82500, .85000, .87500, .90125, .92750 | |
1305 | c *, .96000 | |
1306 | c */ | |
1307 | c data pi4bmx/ | |
1308 | c * .29854, .46181, .54115, .54115, .51709 | |
1309 | c *, .72031, .80976, .96078, .90798, 1.07788 | |
1310 | c *, 1.37389, 1.70101, 2.02716, 1.83687, 1.52644 | |
1311 | c *, 1.57469, 1.35640, .93390, .98531, .84062 | |
1312 | c *, .79188 | |
1313 | c */ | |
1314 | c data pi5ecm/0.75,0.80,0.85,0.90,0.95,1.00,1.05,1.10,1.15/ | |
1315 | c data pi5bmx/0.4,0.5,0.7,0.9,1.1,1.1,0.7,0.6,0.5/ | |
1316 | c | |
1317 | c call utpri('jintcs',ish,ishini,7) | |
1318 | c | |
1319 | c ics=0 | |
1320 | c | |
1321 | c if(idptl(i).gt.10000)then | |
1322 | c call idquad(i,nqi,nai,jci) | |
1323 | c if(nqi.eq.0)then | |
1324 | c idi=999 | |
1325 | c else | |
1326 | c idi=9999 | |
1327 | c endif | |
1328 | c else | |
1329 | c idi=idptl(i) | |
1330 | c endif | |
1331 | c | |
1332 | c if(idptl(j).gt.10000)then | |
1333 | c call idquad(j,nqj,naj,jcj) | |
1334 | c if(nqj.eq.0)then | |
1335 | c idj=999 | |
1336 | c else | |
1337 | c idj=9999 | |
1338 | c endif | |
1339 | c else | |
1340 | c idj=idptl(j) | |
1341 | c endif | |
1342 | c | |
1343 | c do k=1,nflav | |
1344 | c kc(k)=jc(k,1)-jc(k,2) | |
1345 | c if(nq.lt.0)kc(k)=-kc(k) | |
1346 | c enddo | |
1347 | c if(nq.lt.0)nq=-nq | |
1348 | c | |
1349 | c if(ish.ge.7)then | |
1350 | c write(ifch,*)'i:',i,' id_i:',idptl(i),' j:',j,' id_j:',idptl(j) | |
1351 | c write(ifch,*)'E_cm:',ecm,' jc:',(jc(k,1),k=1,6),(jc(k,2),k=1,6) | |
1352 | c write(ifch,*)'nq:',nq,' kc:',(kc(k),k=1,6) | |
1353 | c endif | |
1354 | c | |
1355 | cc check minimal kinetic energy | |
1356 | c ekin=ecm-pptl(5,i)-pptl(5,j) | |
1357 | c if(ekin.lt.amimel)goto1000 | |
1358 | c | |
1359 | cc ------------------------------------------------------------------------- | |
1360 | c if(iabs(idi).gt.1000.or.iabs(idj).gt.1000)then !baryon involved | |
1361 | cc ------------------------------------------------------------------------- | |
1362 | c | |
1363 | c if(nq.eq.6)then !------------baryon-baryon -----> | |
1364 | c | |
1365 | c if(kc(1).eq.kc(2))then !isospin_z=0 | |
1366 | cc pn | |
1367 | c if(ish.ge.7)write(ifch,*)'sig_pn chosen' | |
1368 | c if(ecm.lt.2)then | |
1369 | c call utindx(npn,pnecm,ecm,iecm) | |
1370 | c bmx=pnbmx(iecm) | |
1371 | c else | |
1372 | c p=ecm*sqrt((ecm**2/4./.94**2)-1.) | |
1373 | c if(p.lt.1.)then | |
1374 | c sig=33.+196.*(abs(0.95-p))**2.5 | |
1375 | c elseif(p.lt.2)then | |
1376 | c sig=24.2+8.9*p | |
1377 | c else | |
1378 | c sig=47.3+.513*(alog(p)**2)-4.27*alog(p) | |
1379 | c endif | |
1380 | c bmx=sqrt(sig/10./pi) | |
1381 | c endif | |
1382 | c else !isospin_z=+-1 | |
1383 | cc pp | |
1384 | c if(ish.ge.7)write(ifch,*)'sig_pp chosen' | |
1385 | c if(ecm.le.1.882)then | |
1386 | c bmx=ppbmx(1) | |
1387 | c if(ish.ge.7)write(ifch,*)'b_mx:',bmx | |
1388 | c elseif(ecm.le.1.887)then | |
1389 | c bmx=ppbmx(2) | |
1390 | c if(ish.ge.7)write(ifch,*)'b_mx:',bmx | |
1391 | c elseif(ecm.le.1.893)then | |
1392 | c bmx=ppbmx(3) | |
1393 | c if(ish.ge.7)write(ifch,*)'b_mx:',bmx | |
1394 | c elseif(ecm.le.1.899)then | |
1395 | c bmx=ppbmx(4) | |
1396 | c elseif(ecm.le.1.909)then | |
1397 | c bmx=ppbmx(5) | |
1398 | c elseif(ecm.le.1.913)then | |
1399 | c bmx=ppbmx(6) | |
1400 | c elseif(ecm.le.1.918)then | |
1401 | c bmx=ppbmx(7) | |
1402 | c else | |
1403 | c p=ecm*sqrt((ecm**2/4./.94**2)-1.) | |
1404 | c if(p.lt.0.8)then | |
1405 | c sig=23.5+1000.*(0.7-p)**4 | |
1406 | c elseif(p.lt.1.5)then | |
1407 | c sig=23.5+24.6/(1+exp(-(p-1.2)/0.1)) | |
1408 | c elseif(p.lt.5)then | |
1409 | c sig=41.+60.*(p-0.9)*exp(-1.2*p) | |
1410 | c else | |
1411 | c sig=48+.522*(alog(p)**2)-4.51*alog(p) | |
1412 | c endif | |
1413 | c bmx=sqrt(sig/10./pi) | |
1414 | c endif | |
1415 | c endif | |
1416 | c | |
1417 | c elseif(nq.eq.0)then !-----------baryon-antibaryon ----> | |
1418 | c | |
1419 | c if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0)then !p-ap, n-an | |
1420 | cc app | |
1421 | c if(ish.ge.7)write(ifch,*)'sig_app chosen' | |
1422 | c if(ecm.lt.2)then | |
1423 | c call utindx(napp,appecm,ecm,iecm) | |
1424 | c bmx=appbmx(iecm) | |
1425 | c else | |
1426 | c p=ecm*sqrt((ecm**2/4./.94**2)-1.) | |
1427 | c sig=38.4+77.6*p**(-.64)+.26*(alog(p)**2)-1.2*alog(p) | |
1428 | c bmx=sqrt(sig/10./pi) | |
1429 | c endif | |
1430 | c else !p-an, n-ap | |
1431 | cc apn | |
1432 | c if(ish.ge.7)write(ifch,*)'sig_apn chosen' | |
1433 | c if(ecm.lt.2)then | |
1434 | c call utindx(napn,apnecm,ecm,iecm) | |
1435 | c bmx=apnbmx(iecm) | |
1436 | c else | |
1437 | c p=ecm*sqrt((ecm**2/4./.94**2)-1.) | |
1438 | c sig=133.6*p**(-.7)-1.22*(alog(p)**2)+13.7*alog(p) | |
1439 | c bmx=sqrt(sig/10./pi) | |
1440 | c endif | |
1441 | c endif | |
1442 | c | |
1443 | c elseif(nq.eq.3)then !----------baryon-meson ----> | |
1444 | c | |
1445 | c if(kc(3).eq.0)then | |
1446 | cc no kaons involved(except for K-L) | |
1447 | c if(kc(1).eq.0.or.kc(2).eq.0)then | |
1448 | cc pip | |
1449 | c if(ish.ge.7)write(ifch,*)'sig_pip chosen' | |
1450 | c if(ecm.lt.2.5)then | |
1451 | c call utindx(npip,pipecm,ecm,iecm) | |
1452 | c bmx=pipbmx(iecm) | |
1453 | c else | |
1454 | c p=sqrt(((ecm**2-0.94**2-0.14**2)/2./0.94)**2.-0.14**2) | |
1455 | c sig=16.4+19.3*p**(-.42)+.19*(alog(p)**2) | |
1456 | c bmx=sqrt(sig/10./pi) | |
1457 | c endif | |
1458 | c else | |
1459 | cc pim | |
1460 | c if(ish.ge.7)write(ifch,*)'sig_pim chosen' | |
1461 | c if(ecm.lt.2.5)then | |
1462 | c call utindx(npim,pimecm,ecm,iecm) | |
1463 | c bmx=pimbmx(iecm) | |
1464 | c else | |
1465 | c p=sqrt(((ecm**2-0.94**2-0.14**2)/2./0.94)**2.-0.14**2) | |
1466 | c sig=33.0+14.0*p**(-1.36)+.456*(alog(p)**2)-4.03*alog(p) | |
1467 | c bmx=sqrt(sig/10./pi) | |
1468 | c endif | |
1469 | c endif | |
1470 | c elseif(kc(3).eq.1)then | |
1471 | cc strange particles involved | |
1472 | c if(kc(1).eq.0.or.kc(2).eq.0)then | |
1473 | cc kmn | |
1474 | c if(ish.ge.7)write(ifch,*)'sig_kmn chosen' | |
1475 | c if(ecm.lt.2.5)then | |
1476 | c call utindx(nkmn,kmnecm,ecm,iecm) | |
1477 | c bmx=kmnbmx(iecm) | |
1478 | c else | |
1479 | c p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2) | |
1480 | c sig=25.2+.38*(alog(p)**2)-2.9*alog(p) | |
1481 | c bmx=sqrt(sig/10./pi) | |
1482 | c endif | |
1483 | c else | |
1484 | cc kmp | |
1485 | c if(ish.ge.7)write(ifch,*)'sig_kmp chosen' | |
1486 | c if(ecm.lt.2.5)then | |
1487 | c call utindx(nkmp,kmpecm,ecm,iecm) | |
1488 | c bmx=kmpbmx(iecm) | |
1489 | c else | |
1490 | c p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2) | |
1491 | c sig=32.1+.66*(alog(p)**2)-5.6*alog(p) | |
1492 | c bmx=sqrt(sig/10./pi) | |
1493 | c endif | |
1494 | c endif | |
1495 | c elseif(kc(3).eq.-1)then | |
1496 | cc strange particles involved | |
1497 | c if(kc(1).eq.3.or.kc(2).eq.3)then | |
1498 | cc kpp | |
1499 | c if(ish.ge.7)write(ifch,*)'sig_kpp chosen' | |
1500 | c if(ecm.lt.2.5)then | |
1501 | c call utindx(nkpp,kppecm,ecm,iecm) | |
1502 | c bmx=kppbmx(iecm) | |
1503 | c else | |
1504 | c p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2) | |
1505 | c sig=18.1+.26*(alog(p)**2)-alog(p) | |
1506 | c bmx=sqrt(sig/10./pi) | |
1507 | c endif | |
1508 | c else | |
1509 | cc kpn | |
1510 | c if(ish.ge.7)write(ifch,*)'sig_kpn chosen' | |
1511 | c if(ecm.lt.2.5)then | |
1512 | c call utindx(nkpn,kpnecm,ecm,iecm) | |
1513 | c bmx=kpnbmx(iecm) | |
1514 | c else | |
1515 | c p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2) | |
1516 | c sig=18.7+.21*(alog(p)**2)-.89*alog(p) | |
1517 | c bmx=sqrt(sig/10./pi) | |
1518 | c endif | |
1519 | c endif | |
1520 | c elseif(kc(3).ge.2)then | |
1521 | cc two strange particles involved | |
1522 | c bmx=1. | |
1523 | c elseif(kc(3).le.-2)then | |
1524 | cc two strange particles involved | |
1525 | c bmx=1. | |
1526 | c endif | |
1527 | c | |
1528 | c else !-----------baryon_cluster-anything----> | |
1529 | c | |
1530 | c write(ifch,*)'i:',i,' id_i:',idptl(i),' j:',j,' id_j:',idptl(j) | |
1531 | c write(ifch,*)'r_i:',radptl(i),' r_j:',radptl(j) | |
1532 | c write(ifch,*)'E_cm:',ecm,' jc:',(jc(k,1),k=1,6),(jc(k,2),k=1,6) | |
1533 | c write(ifch,*)'nq:',nq,' kc:',(kc(k),k=1,6) | |
1534 | c bmx=radptl(i)+radptl(j) | |
1535 | c | |
1536 | c endif ! <-------------- | |
1537 | c | |
1538 | cc ------------------------------------- | |
1539 | c else ! meson-meson | |
1540 | cc ------------------------------------- | |
1541 | c | |
1542 | c call idquad(i,nqi,nai,jci) | |
1543 | c call idquad(j,nqj,naj,jcj) | |
1544 | c do l=1,nflav | |
1545 | c kci(l)=jci(l,1)-jci(l,2) | |
1546 | c kcj(l)=jcj(l,1)-jcj(l,2) | |
1547 | c enddo | |
1548 | c | |
1549 | c if(kci(3).eq.0.and.pptl(5,i).le.1.0 | |
1550 | c * .and.kcj(3).eq.0.and.pptl(5,j).le.1.0)then | |
1551 | c | |
1552 | c if(kci(1).eq.0.and.kci(2).eq.0.and.pptl(5,i).le.0.5)then | |
1553 | c if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then | |
1554 | cc pi0 pi0 | |
1555 | c goto104 | |
1556 | c elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then | |
1557 | cc pi0 pi+ | |
1558 | c goto102 | |
1559 | c elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then | |
1560 | cc pi0 pi- | |
1561 | c goto103 | |
1562 | c else | |
1563 | cc pi0 eta | |
1564 | c goto105 | |
1565 | c endif | |
1566 | c elseif(kci(1).eq.1.and.kci(2).eq.-1)then | |
1567 | c if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then | |
1568 | cc pi+ pi0 | |
1569 | c goto102 | |
1570 | c elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then | |
1571 | cc pi+ pi+ | |
1572 | c goto101 | |
1573 | c elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then | |
1574 | cc pi+ pi- | |
1575 | c goto104 | |
1576 | c else | |
1577 | cc pi+ eta | |
1578 | c goto105 | |
1579 | c endif | |
1580 | c elseif(kci(1).eq.-1.and.kci(2).eq.1)then | |
1581 | c if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then | |
1582 | cc pi- pi0 | |
1583 | c goto103 | |
1584 | c elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then | |
1585 | cc pi- pi+ | |
1586 | c goto104 | |
1587 | c elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then | |
1588 | cc pi- pi- | |
1589 | c goto101 | |
1590 | c else | |
1591 | cc pi- eta | |
1592 | c goto105 | |
1593 | c endif | |
1594 | c else | |
1595 | c if(pptl(5,j).le.0.5)then | |
1596 | cc eta pi | |
1597 | c goto105 | |
1598 | c else | |
1599 | cc eta eta | |
1600 | c goto106 | |
1601 | c endif | |
1602 | c endif | |
1603 | c101 continue | |
1604 | c if(ish.ge.7)write(ifch,*)'sig_pi+pi+ chosen' | |
1605 | c call utindx(npi1,pi1ecm,ecm,iecm) | |
1606 | c bmx=pi1bmx(iecm) | |
1607 | c goto110 | |
1608 | c102 continue | |
1609 | c if(ish.ge.7)write(ifch,*)'sig_pi+pi0 chosen' | |
1610 | c call utindx(npi2,pi2ecm,ecm,iecm) | |
1611 | c bmx=pi2bmx(iecm) | |
1612 | c goto110 | |
1613 | c103 continue | |
1614 | c if(ish.ge.7)write(ifch,*)'sig_pi-pi0 chosen' | |
1615 | c call utindx(npi3,pi3ecm,ecm,iecm) | |
1616 | c bmx=pi3bmx(iecm) | |
1617 | c goto110 | |
1618 | c104 continue | |
1619 | c if(ish.ge.7)write(ifch,*)'sig_pi-pi+ chosen' | |
1620 | c call utindx(npi4,pi4ecm,ecm,iecm) | |
1621 | c bmx=pi4bmx(iecm) | |
1622 | c goto110 | |
1623 | c105 continue | |
1624 | c if(ish.ge.7)write(ifch,*)'sig_pi_eta chosen' | |
1625 | c call utindx(npi5,pi5ecm,ecm,iecm) | |
1626 | c bmx=pi5bmx(iecm) | |
1627 | c goto110 | |
1628 | c106 continue | |
1629 | c if(ish.ge.7)write(ifch,*)'sig_eta_eta chosen' | |
1630 | c bmx=0.4 ! approx. sig=5mb | |
1631 | c110 continue | |
1632 | c | |
1633 | c else !something else involved, strange etc. | |
1634 | c | |
1635 | c bmx=0.6 ! approx. sig=10mb | |
1636 | c | |
1637 | c endif | |
1638 | c | |
1639 | cc -------------------------------- | |
1640 | c endif | |
1641 | cc -------------------------------- | |
1642 | c | |
1643 | c if(bij.le.bmx)ics=1 | |
1644 | c if(ish.ge.7)then | |
1645 | c write(ifch,*)'b_mx:',bmx,' b_ij:',bij,' ics:',ics | |
1646 | c endif | |
1647 | c | |
1648 | c1000 continue | |
1649 | c call utprix('jintcs',ish,ishini,7) | |
1650 | c return | |
1651 | c end | |
1652 | c | |
1653 | cc---------------------------------------------------------------------- | |
1654 | c subroutine jintel(i,j,p,amim,xaver) | |
1655 | cc---------------------------------------------------------------------- | |
1656 | cc elastic scattering of ptls i,j | |
1657 | cc---------------------------------------------------------------------- | |
1658 | c include 'epos.inc' | |
1659 | c real xaver(4) | |
1660 | c real p(5),u(5),pei(5),pej(5) | |
1661 | c | |
1662 | cc determine momenta of outgoing particles (pei,pej) | |
1663 | cc ------------------------------------------------- | |
1664 | c if(p(5).le.(pptl(5,i)+pptl(5,j))*.99)then | |
1665 | c if(ish.ge.1)then | |
1666 | c call utmsg('jintel') | |
1667 | c write(ifch,132)p(5),pptl(5,i)+pptl(5,j) | |
1668 | c132 format(1x,'***** m_fus < m_i+m_j ---> qcm set zero ( ' | |
1669 | c *,2f10.3,' )') | |
1670 | c write(ifch,133)'p_i: ',(pptl(k,i),k=1,5) | |
1671 | c write(ifch,133)'p_j: ',(pptl(k,j),k=1,5) | |
1672 | c write(ifch,133)'p_fus:',(p(k),k=1,5) | |
1673 | c133 format(1x,a6,3x,5f10.3) | |
1674 | c call utmsgf | |
1675 | c endif | |
1676 | c qcm=0. | |
1677 | c elseif(p(5).le.pptl(5,i)+pptl(5,j))then | |
1678 | c qcm=0. | |
1679 | c else | |
1680 | c qcm=utpcm(p(5),pptl(5,i),pptl(5,j)) | |
1681 | c endif | |
1682 | c | |
1683 | cc isotropic | |
1684 | c u(3)=2.*rangen()-1. | |
1685 | c phi=2.*pi*rangen() | |
1686 | c u(1)=sqrt(1.-u(3)**2)*cos(phi) | |
1687 | c u(2)=sqrt(1.-u(3)**2)*sin(phi) | |
1688 | c do 47 k=1,3 | |
1689 | c pei(k)= qcm*u(k) | |
1690 | c47 pej(k)=-qcm*u(k) | |
1691 | c | |
1692 | cc nonisotropic | |
1693 | cc-c pt=sqrt(2./pi)*ptq*sqrt(-2*alog(rangen()) !2-dim Gauss | |
1694 | cc-c if(pt.ge.qcm)pt=rangen()*qcm | |
1695 | cc-c qpl=sqrt(qcm**2-pt**2) | |
1696 | cc-c u(3)=qpl | |
1697 | cc-c phi=2.*pi*rangen() | |
1698 | cc-c u(1)=pt*cos(phi) | |
1699 | cc-c u(2)=pt*sin(phi) | |
1700 | cc-c call utaxis(i,j,a1,a2,a3) | |
1701 | cc-c ivt=1 | |
1702 | cc-c if(a3.lt.0.)then | |
1703 | cc-c a1=-a1 | |
1704 | cc-c a2=-a2 | |
1705 | cc-c a3=-a3 | |
1706 | cc-c ivt=-1 | |
1707 | cc-c endif | |
1708 | cc-c call utrota(-1,a1,a2,a3,u(1),u(2),u(3)) | |
1709 | cc-c do 47 k=1,3 | |
1710 | cc-c pei(k)= u(k)*ivt | |
1711 | cc-c47 pej(k)=-u(k)*ivt | |
1712 | c | |
1713 | c pei(4)=sqrt(qcm**2+pptl(5,i)**2) | |
1714 | c pej(4)=sqrt(qcm**2+pptl(5,j)**2) | |
1715 | c pei(5)=pptl(5,i) | |
1716 | c pej(5)=pptl(5,j) | |
1717 | c call utlobo(-1,p(1),p(2),p(3),p(4),p(5) | |
1718 | c *,pei(1),pei(2),pei(3),pei(4)) | |
1719 | c call utlobo(-1,p(1),p(2),p(3),p(4),p(5) | |
1720 | c *,pej(1),pej(2),pej(3),pej(4)) | |
1721 | c | |
1722 | cc fill /cptl/ | |
1723 | cc ----------- | |
1724 | c do 49 lo=1,2 | |
1725 | c nptl=nptl+1 | |
1726 | c if(lo.eq.1)ij=i | |
1727 | c if(lo.eq.2)ij=j | |
1728 | c do 48 k=1,5 | |
1729 | c if(lo.eq.1)pptl(k,nptl)=pei(k) | |
1730 | c if(lo.eq.2)pptl(k,nptl)=pej(k) | |
1731 | c48 continue | |
1732 | c istptl(nptl)=0 | |
1733 | c idptl(nptl)=idptl(ij) | |
1734 | c ibptl(1,nptl)=ibptl(1,ij) | |
1735 | c ibptl(2,nptl)=ibptl(2,ij) | |
1736 | c ibptl(3,nptl)=ibptl(3,ij) | |
1737 | c ibptl(4,nptl)=ibptl(4,ij) | |
1738 | c xorptl(1,nptl)=xaver(1) | |
1739 | c xorptl(2,nptl)=xaver(2) | |
1740 | c xorptl(3,nptl)=xaver(3) | |
1741 | c xorptl(4,nptl)=xaver(4) | |
1742 | c iorptl(nptl)=i | |
1743 | c jorptl(nptl)=j | |
1744 | c tivptl(1,nptl)=xaver(4) | |
1745 | c call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm) | |
1746 | c tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen())) | |
1747 | c ifrptl(1,nptl)=0 | |
1748 | c ifrptl(2,nptl)=0 | |
1749 | c ityptl(nptl)=ityptl(ij) | |
1750 | c49 continue | |
1751 | c | |
1752 | c1000 return | |
1753 | c end | |
1754 | c | |
1755 | cc---------------------------------------------------------------------- | |
1756 | c subroutine jintfs(i,j,p,nq,jc,amim,iret) | |
1757 | cc---------------------------------------------------------------------- | |
1758 | cc input: | |
1759 | cc i,j: particle indices | |
1760 | cc output: | |
1761 | cc p: 5-momentum of fused ptl | |
1762 | cc nq: net quark number of fused ptl | |
1763 | cc jc: jc of fused ptl | |
1764 | cc amim : minimum mass of fused ptl | |
1765 | cc iret: return code | |
1766 | cc 0=ok | |
1767 | cc 1=mass p(5) less than amim | |
1768 | cc---------------------------------------------------------------------- | |
1769 | c include 'epos.inc' | |
1770 | c integer jci(nflav,2),jcj(nflav,2),jc(nflav,2) | |
1771 | c real p(5) | |
1772 | c double precision ppfus(4),pp52 | |
1773 | c | |
1774 | c iret=0 | |
1775 | c | |
1776 | c do 35 k=1,4 | |
1777 | c p(k)=pptl(k,i)+pptl(k,j) | |
1778 | c35 ppfus(k)=dble(p(k)) | |
1779 | c pp52=ppfus(4)**2-ppfus(3)**2-ppfus(2)**2-ppfus(1)**2 | |
1780 | c if(pp52.le.0.)then | |
1781 | c if(ish.ge.1)then | |
1782 | c call utmsg('jintfs') | |
1783 | c write(ifch,*)'***** mfus**2 < 0 (',pp52,' )' | |
1784 | c write(ifch,*)(ppfus(m),m=1,4) | |
1785 | c call utmsgf | |
1786 | c endif | |
1787 | c goto1001 | |
1788 | c endif | |
1789 | c p(5)=sngl(dsqrt(pp52)) | |
1790 | c | |
1791 | c call idquad(i,nqi,nai,jci) | |
1792 | c call idquad(j,nqj,naj,jcj) | |
1793 | c do 29 n=1,nflav | |
1794 | c do 29 k=1,2 | |
1795 | c29 jc(n,k)=jci(n,k)+jcj(n,k) | |
1796 | c nq=0 | |
1797 | c do 54 n=1,nflav | |
1798 | c54 nq=nq+jc(n,1)-jc(n,2) | |
1799 | c | |
1800 | c call idcomj(jc) | |
1801 | c amim=utamnz(jc,5)+amimfs | |
1802 | c if(p(5).lt.amim)goto1001 | |
1803 | c goto1000 | |
1804 | c1001 iret=1 | |
1805 | c1000 return | |
1806 | c end | |
1807 | c | |
1808 | cc---------------------------------------------------------------------- | |
1809 | c subroutine jintfu(i,j,p,jc) | |
1810 | cc---------------------------------------------------------------------- | |
1811 | cc input: | |
1812 | cc i,j: particle indices | |
1813 | cc p,jc: momentum and jc of fused object | |
1814 | cc outout: | |
1815 | cc id, ib, ist, ior, jor, ifr, ity of fused particle | |
1816 | cc written to /cptl/ after increasing nptl by 1 | |
1817 | cc---------------------------------------------------------------------- | |
1818 | c | |
1819 | c include 'epos.inc' | |
1820 | c parameter (mxlook=10000,mxdky=2000) | |
1821 | c common/dkytab/look(mxlook),cbr(mxdky),mode(5,mxdky) | |
1822 | c real p(5) | |
1823 | c integer jc(nflav,2),ic(2),ib(4) | |
1824 | c | |
1825 | c nptl=nptl+1 | |
1826 | c | |
1827 | cc momentum | |
1828 | c do k=1,5 | |
1829 | c pptl(k,nptl)=p(k) | |
1830 | c enddo | |
1831 | c amf=p(5) | |
1832 | c | |
1833 | cc determine idr, ib(1-4) | |
1834 | c idr=0 | |
1835 | c do 40 nf=1,nflav | |
1836 | c do 40 ij=1,2 | |
1837 | c if(jc(nf,ij).ge.10)idr=7*10**8 | |
1838 | c40 continue | |
1839 | c if(idr/10**8.ne.7)then | |
1840 | c call idenco(jc,ic,ireten) | |
1841 | c if(ireten.eq.1)call utstop('jintfu: idenco ret code = 1&') | |
1842 | c id=idtra(ic,0,0,3) | |
1843 | c43 amc=amf | |
1844 | c call idres(id,amc,idr,iadj,1) | |
1845 | c if(idr.ne.0)then | |
1846 | c lid=look(iabs(idr)) | |
1847 | c if((lid.le.0.or.lid.gt.0.and.mode(2,lid).eq.0) | |
1848 | c *.and.pptl(5,nptl).gt.amc+1e-3)then | |
1849 | c amf=amf+0.010 | |
1850 | c goto43 | |
1851 | c endif | |
1852 | c if((lid.le.0.or.lid.gt.0.and.mode(2,lid).eq.0) | |
1853 | c *.and.abs(amc-pptl(5,nptl)).gt.1e-3)then | |
1854 | c if(ish.ge.1)then | |
1855 | c call utmsg('jintfu') | |
1856 | c write(ifch,*)'***** not on mass shell after fusion: ' | |
1857 | c *,pptl(5,nptl),amc | |
1858 | c call utmsgf | |
1859 | c endif | |
1860 | c endif | |
1861 | c endif | |
1862 | c if(idr.eq.0)then | |
1863 | c if(mod(ic(1),100).ne.0.or.mod(ic(2),100).ne.0)then | |
1864 | c idr=9*10**8 | |
1865 | c else | |
1866 | c idr=8*10**8+ic(1)*100+ic(2)/100 | |
1867 | c endif | |
1868 | c endif | |
1869 | c else | |
1870 | c call idtrbi(jc,ib(1),ib(2),ib(3),ib(4)) | |
1871 | c idr=idr | |
1872 | c *+mod(jc(1,1)+jc(2,1)+jc(3,1)+jc(4,1),10**4)*10**4 | |
1873 | c *+mod(jc(1,2)+jc(2,2)+jc(3,2)+jc(4,2),10**4) | |
1874 | c if(ish.ge.7)write(ifch,*)'ib:',(ib(kk),kk=1,4) | |
1875 | c endif | |
1876 | c | |
1877 | c n=nptl | |
1878 | c | |
1879 | cc fill /cptl/ | |
1880 | cc ----------- | |
1881 | c idptl(n)=idr | |
1882 | c do k=1,4 | |
1883 | c ibptl(k,n)=ib(k) | |
1884 | c enddo | |
1885 | c istptl(n)=0 | |
1886 | c iorptl(n)=i | |
1887 | c jorptl(n)=j | |
1888 | c ifrptl(1,n)=0 | |
1889 | c ifrptl(2,n)=0 | |
1890 | c ityptl(n)=50 | |
1891 | c | |
1892 | c return | |
1893 | c end | |
1894 | c | |
1895 | c---------------------------------------------------------------------- | |
1896 | subroutine jintpo | |
1897 | c---------------------------------------------------------------------- | |
1898 | c parton-ladder-fusion -- 3D grid -- based on string segments | |
1899 | c---------------------------------------------------------------------- | |
1900 | include 'epos.inc' | |
1901 | include 'epos.incico' | |
1902 | common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl) | |
1903 | *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3) | |
1904 | parameter(xgrid=8,mxgrid=10,kgrid=6,kegrid=3) | |
1905 | parameter(ssgrid=0.5*xgrid*kgrid*kegrid) | |
1906 | parameter(mssgrid=mxgrid*kgrid*kegrid) | |
1907 | parameter(mxcl=1000,mxcli=200) | |
1908 | integer idropgrid(mxgrid,mxgrid,mssgrid) | |
1909 | & ,jdropgrid(mxgrid,mxgrid,mssgrid) | |
1910 | & ,jclu(mssgrid),nsegp4(mxcl) | |
1911 | & ,irep(mxcl),mmji(mxcl,mxcli) | |
1912 | & ,jccl(mxcl,nflav,2),nseg(mxcl),mseg(mxcl),kclu(mxcl) | |
1913 | & ,naseg(0:mxcl),nfseg(mxcl),nsegmx(mxcl) | |
1914 | & ,ic(2),jc(nflav,2),ke(6),jcjj(nflav,2) | |
1915 | & ,nsegmt(mxptl) !,nclk(mssgrid) | |
1916 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
1917 | & ,pptld(5,mxcl),p4tmp | |
1918 | common /cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat | |
1919 | common/cdelzet/delzet,delsgr /cvocell/vocell | |
1920 | common/cranphi/ranphi,ranecc,weiecc | |
1921 | logical first | |
1922 | double precision ptest(4),ttest,p52,p4mean(mxcl) !,ppp(5),qqq(5),www(5) | |
1923 | save ncntpo | |
1924 | data ncntpo /0/ | |
1925 | ncntpo=ncntpo+1 | |
1926 | call utpri('jintpo',ish,ishini,4) | |
1927 | ||
1928 | if(ttaus.gt.kgrid)then | |
1929 | write(ifmt,*)'ttaus,kgrid :',ttaus,kgrid | |
1930 | stop'jintpo: ttaus too big.' | |
1931 | endif | |
1932 | tauico=ttaus | |
1933 | fegrid=yhaha/5.36 !rapidity range/rap range at RHIC | |
1934 | sgrid=ttaus*ssgrid/kgrid *fegrid | |
1935 | msgrid=max(1.4d0,ttaus)*mssgrid/kgrid *fegrid | |
1936 | if(mod(msgrid,2).eq.0)msgrid=msgrid+1 | |
1937 | vocell= 2*xgrid/mxgrid * 2*xgrid/mxgrid * 2*sgrid/msgrid | |
1938 | delzet=2*sgrid/msgrid/ttaus | |
1939 | delsgr=2*sgrid/msgrid | |
1940 | nptla=nptl | |
1941 | nsegsuj=max(3,nint(float(nsegsu)*fegrid)) | |
1942 | !print*,'dx,dz:',2*xgrid/mxgrid,2*sgrid/msgrid | |
1943 | !. ,' nx nz:', mxgrid,msgrid, ttaus | |
1944 | ||
1945 | c...compute x,y,z | |
1946 | ||
1947 | if(ish.ge.6)write(ifch,*)'compute x,y,z' | |
1948 | do n=1,nptla | |
1949 | if(istptl(n).eq.0)then | |
1950 | call jtain(n,xptl(n),yptl(n),zptl(n),tptl(n),nnn,0) | |
1951 | call jtaus(zptl(n),tz,sptl(n)) | |
1952 | strap=1e10 | |
1953 | xpl=tptl(n)+zptl(n) | |
1954 | xmi=tptl(n)-zptl(n) | |
1955 | if(xmi.gt.0.0.and.xpl.gt.0.0)strap=0.5*log(xpl/xmi) | |
1956 | dezptl(n)=strap !space-time-rapidity | |
1957 | c if(iorptl(n).eq.4142)write(ifch,*)'ini',n | |
1958 | c & ,pptl(4,n),pptl(3,n),idptl(n),xptl(n),yptl(n),zptl(n),sptl(n) | |
1959 | else | |
1960 | xptl(n)=0. | |
1961 | yptl(n)=0. | |
1962 | zptl(n)=0. | |
1963 | sptl(n)=0. | |
1964 | dezptl(n)=1e10 | |
1965 | endif | |
1966 | enddo | |
1967 | ntry=0 | |
1968 | ||
1969 | c...valid particles | |
1970 | ||
1971 | do n=1,nptla | |
1972 | iaaptl(n)=1 | |
1973 | pt2=pptl(1,n)**2+pptl(2,n)**2 | |
1974 | am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2 | |
1975 | if(istptl(n).eq.0)then | |
1976 | call idquac(n,idum1,idum2,idum3,jc) | |
1977 | else | |
1978 | jc(4,1)=0 | |
1979 | jc(4,2)=0 | |
1980 | endif | |
1981 | if(istptl(n).ne.0.or.ityptl(n).ge.60)then | |
1982 | iaaptl(n)=0 | |
1983 | !elseif(ityptl(n).eq.41.or.ityptl(n).eq.51)then | |
1984 | ! iaaptl(n)=0 !avoid particles coming from remnant droplet decay | |
1985 | ! !(already droplet decay products !) | |
1986 | elseif(abs(am2tmp-pptl(5,n)*pptl(5,n)).gt.30.)then | |
1987 | iaaptl(n)=0 !to discard off shell particles | |
1988 | elseif(jc(4,1).ne.0.or.jc(4,2).ne.0)then | |
1989 | iaaptl(n)=0 !to discard charmed particles | |
1990 | !elseif(pt2.lt.1.e-3)then | |
1991 | ! iaaptl(n)=0 !to discard slow proton (spectators) | |
1992 | elseif(pt2.gt.(1.5*ptclu)**2)then | |
1993 | if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)iaaptl(n)=0 | |
1994 | elseif(pt2.gt.(0.5*ptclu)**2)then | |
1995 | if(rangen().lt.(sqrt(pt2)-0.5*ptclu)/ptclu)then | |
1996 | if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)iaaptl(n)=0 | |
1997 | endif | |
1998 | endif | |
1999 | enddo | |
2000 | ||
2001 | c...Start cluster formation | |
2002 | ||
2003 | 8888 continue | |
2004 | ntry=ntry+1 | |
2005 | if(ntry.gt.90) !do not put more than 100 or change limit for p4mean | |
2006 | &call utstop('jintpo: cluster formation impossible ! &') | |
2007 | nptl=nptla | |
2008 | ||
2009 | do 1 k=1,msgrid | |
2010 | do 1 j=1,mxgrid | |
2011 | do 1 i=1,mxgrid | |
2012 | idropgrid(i,j,k)=0 | |
2013 | 1 continue | |
2014 | ||
2015 | c...count string segments in cell | |
2016 | ||
2017 | if(ish.ge.6)write(ifch,*)'count string segments in cell' | |
2018 | do n=1,nptla | |
2019 | if(iaaptl(n).ne.0)then | |
2020 | i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2021 | j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2022 | k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid) | |
2023 | if( i.ge.1.and.i.le.mxgrid | |
2024 | & .and.j.ge.1.and.j.le.mxgrid | |
2025 | & .and.k.ge.1.and.k.le.msgrid) | |
2026 | & idropgrid(i,j,k)=idropgrid(i,j,k)+1 | |
2027 | endif | |
2028 | enddo | |
2029 | ||
2030 | c...check high pt segments | |
2031 | c !to use this part one has to define: | |
2032 | c !... 1 = valid | |
2033 | c !... -1 = valid but high pt | |
2034 | c !... 0 = not valid | |
2035 | c | |
2036 | c if(ish.ge.6)write(ifch,*)'check high pt segments' | |
2037 | c do n=1,nptla | |
2038 | c if(iaaptl(n).eq.-1)then | |
2039 | c i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2040 | c j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2041 | c k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid) | |
2042 | c if( i.ge.1.and.i.le.mxgrid | |
2043 | c & .and.j.ge.1.and.j.le.mxgrid | |
2044 | c & .and.k.ge.1.and.k.le.msgrid)then | |
2045 | c if(idropgrid(i,j,k).lt.2*nsegce)then !surface segments | |
2046 | c iaaptl(n)=0 | |
2047 | c idropgrid(i,j,k)=idropgrid(i,j,k)-1 | |
2048 | c if(idropgrid(i,j,k).lt.0)stop'jintpo: not possible. ' | |
2049 | c else | |
2050 | c iaaptl(n)=1 | |
2051 | c endif | |
2052 | c endif | |
2053 | c endif | |
2054 | c enddo | |
2055 | ||
2056 | c...identify clusters | |
2057 | ||
2058 | if(ish.ge.6)write(ifch,*)'identify clusters' | |
2059 | do k=1,msgrid !~~~~~~k-loop | |
2060 | jjj=0 | |
2061 | do j=1,mxgrid | |
2062 | first=.true. | |
2063 | do i=1,mxgrid | |
2064 | if(idropgrid(i,j,k).ge.nsegce)then | |
2065 | if(first)then | |
2066 | ifirst=i | |
2067 | jjj=jjj+1 | |
2068 | if(jjj.gt.mxcl)stop'jintpo: mxcl too small. ' | |
2069 | irep(jjj)=0 | |
2070 | jj=jjj | |
2071 | first=.false. | |
2072 | endif | |
2073 | jdropgrid(i,j,k)=jj | |
2074 | if(j.gt.1.and.jdropgrid(i,j-1,k).ne.0)then | |
2075 | jjo=jdropgrid(i,j-1,k) | |
2076 | if(jjo.lt.jj)then | |
2077 | if(jj.eq.jjj)jjj=jjj-1 | |
2078 | jjx=jj | |
2079 | jj=jjo | |
2080 | do ii=ifirst,i | |
2081 | jdropgrid(ii,j,k)=jj | |
2082 | if(jdropgrid(ii,j-1,k).eq.jjx)then | |
2083 | if(jjx.gt.jjj)jjj=jjj+1 | |
2084 | jja=jjx | |
2085 | jjb=jj | |
2086 | 90 continue | |
2087 | if(irep(jja).eq.0.or.irep(jja).eq.jjb)then | |
2088 | irep(jja)=jjb | |
2089 | else | |
2090 | mn=min(irep(jja),jjb) | |
2091 | mx=max(irep(jja),jjb) | |
2092 | irep(jja)=mn | |
2093 | jja=mx | |
2094 | jjb=mn | |
2095 | goto90 | |
2096 | endif | |
2097 | endif | |
2098 | enddo | |
2099 | elseif(jdropgrid(i,j-1,k).gt.jj)then | |
2100 | irep(jjo)=jj | |
2101 | endif | |
2102 | endif | |
2103 | else | |
2104 | jdropgrid(i,j,k)=0 | |
2105 | first=.true. | |
2106 | endif | |
2107 | enddo | |
2108 | enddo | |
2109 | !~~~~cluster jj ---> cluster irep(jj) | |
2110 | do jj=jjj,1,-1 | |
2111 | if(irep(jj).ne.0)then | |
2112 | do j=1,mxgrid | |
2113 | do i=1,mxgrid | |
2114 | if(jdropgrid(i,j,k).eq.jj)jdropgrid(i,j,k)=irep(jj) | |
2115 | enddo | |
2116 | enddo | |
2117 | endif | |
2118 | enddo | |
2119 | !~~~~~remove empty cluster indices | |
2120 | jjjx=jjj | |
2121 | jjj=0 | |
2122 | jj=0 | |
2123 | do jjx=1,jjjx | |
2124 | if(irep(jjx).eq.0)then | |
2125 | jj=jj+1 | |
2126 | jjj=jjj+1 | |
2127 | else | |
2128 | do j=1,mxgrid | |
2129 | do i=1,mxgrid | |
2130 | if(jdropgrid(i,j,k).gt.jj) | |
2131 | & jdropgrid(i,j,k)=jdropgrid(i,j,k)-1 | |
2132 | enddo | |
2133 | enddo | |
2134 | endif | |
2135 | enddo | |
2136 | !~~~~~ | |
2137 | jclu(k)=jjj | |
2138 | enddo !~~~~~~~~~~~~~~~~~ END k-loop | |
2139 | ||
2140 | c...absolute clusters numbering (for all k) | |
2141 | ||
2142 | if(ish.ge.6)write(ifch,*)'absolute clusters numbering' | |
2143 | jjj=jclu(1) | |
2144 | do k=2,msgrid | |
2145 | do j=1,mxgrid | |
2146 | do i=1,mxgrid | |
2147 | if(jdropgrid(i,j,k).gt.0) | |
2148 | & jdropgrid(i,j,k)=jdropgrid(i,j,k)+jjj | |
2149 | enddo | |
2150 | enddo | |
2151 | jjj=jjj+jclu(k) | |
2152 | enddo | |
2153 | jjs=0 | |
2154 | do k=1,msgrid | |
2155 | do jj=1,jclu(k) | |
2156 | jjs=jjs+1 | |
2157 | kclu(jjs)=k | |
2158 | nseg(jjs)=0 | |
2159 | nsegp4(jjs)=0 | |
2160 | p4mean(jjs)=0d0 | |
2161 | enddo | |
2162 | enddo | |
2163 | nptlb=nptl+jjj | |
2164 | ||
2165 | c...calculate mean energy of segments going into in clusters | |
2166 | ||
2167 | if(ish.ge.6)write(ifch,*) | |
2168 | &'calculate mean energy of segments going into in clusters' | |
2169 | do 95 n=1,nptla | |
2170 | if(iaaptl(n).eq.0)goto 95 | |
2171 | i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2172 | j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2173 | k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid) | |
2174 | if( i.ge.1.and.i.le.mxgrid | |
2175 | & .and.j.ge.1.and.j.le.mxgrid | |
2176 | & .and.k.ge.1.and.k.le.msgrid)then | |
2177 | jj=jdropgrid(i,j,k) | |
2178 | if(jj.gt.0)then | |
2179 | nsegp4(jj)=nsegp4(jj)+1 | |
2180 | p4mean(jj)=p4mean(jj)+dble(pptl(4,n)) | |
2181 | c if(jj.eq.9)write(ifch,*)'after',n | |
2182 | c & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k | |
2183 | endif | |
2184 | endif | |
2185 | 95 continue | |
2186 | ||
2187 | do jj=1,jjj | |
2188 | if(nsegp4(jj).gt.5)then | |
2189 | p4mean(jj)=p4mean(jj)/dble(nsegp4(jj)) | |
2190 | else | |
2191 | p4mean(jj)=1d50 | |
2192 | endif | |
2193 | enddo | |
2194 | ||
2195 | c...mark segments going into in clusters, count them | |
2196 | ||
2197 | if(ish.ge.6)write(ifch,*)'mark segments going into in clusters' | |
2198 | do 96 n=1,nptla | |
2199 | if(iaaptl(n).eq.0)goto96 | |
2200 | i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2201 | j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2202 | k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid) | |
2203 | if( i.ge.1.and.i.le.mxgrid | |
2204 | & .and.j.ge.1.and.j.le.mxgrid | |
2205 | & .and.k.ge.1.and.k.le.msgrid)then | |
2206 | jj=jdropgrid(i,j,k) | |
2207 | if(jj.gt.0)then | |
2208 | ! count only particles with reasonnable energy | |
2209 | pt2=pptl(1,n)**2+pptl(2,n)**2+pptl(5,n)**2 | |
2210 | am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2 | |
2211 | if(abs(am2tmp).lt.0.1.or. | |
2212 | & dble(pptl(4,n)).le.100.d0/dble(ntry)*p4mean(jj))then | |
2213 | istptl(n) = 3 | |
2214 | nseg(jj)=nseg(jj)+1 | |
2215 | c if(jj.eq.9)write(ifch,*)'after',n | |
2216 | c & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k | |
2217 | else | |
2218 | c write(ifmt,*)'Rejected particle : ',n,idptl(n) | |
2219 | if(ish.ge.1)write(ifch,*)'Rejected particle : ',n,idptl(n) | |
2220 | & ,ityptl(n),(pptl(k,n),k=1,5),am2tmp | |
2221 | & ,dble(pptl(4,n))/dble(nsegp4(jj)),p4mean(jj),nsegp4(jj),jj | |
2222 | nsegp4(jj)=nsegp4(jj)-1 | |
2223 | iaaptl(n)=0 | |
2224 | if(nsegp4(jj).gt.5)then | |
2225 | p4mean(jj)=(p4mean(jj)*dble(nsegp4(jj)+1)-pptl(4,n)) | |
2226 | & /dble(nsegp4(jj)) | |
2227 | else | |
2228 | p4mean(jj)=1d50 | |
2229 | endif | |
2230 | endif | |
2231 | endif | |
2232 | endif | |
2233 | 96 continue | |
2234 | ||
2235 | ||
2236 | c...add segments moving towards clusters | |
2237 | ||
2238 | if(ish.ge.6)write(ifch,*)'add segments moving towards clusters' | |
2239 | if(iocluin.eq.1)then | |
2240 | do 93 n=1,nptla | |
2241 | if(iaaptl(n).eq.0)goto93 | |
2242 | ihit=0 | |
2243 | i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2244 | j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2245 | k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid) | |
2246 | if( i.ge.1.and.i.le.mxgrid | |
2247 | & .and.j.ge.1.and.j.le.mxgrid | |
2248 | & .and.k.ge.1.and.k.le.msgrid)then | |
2249 | jgr=jdropgrid(i,j,k) | |
2250 | if(jgr.eq.0)then | |
2251 | if(i.ge.mxgrid/2)then | |
2252 | isi=-1 | |
2253 | else | |
2254 | isi=1 | |
2255 | endif | |
2256 | if(j.ge.mxgrid/2)then | |
2257 | jsi=-1 | |
2258 | else | |
2259 | jsi=1 | |
2260 | endif | |
2261 | do ii=i,i+2*isi,isi | |
2262 | do jj=j,j+2*jsi,jsi | |
2263 | if(.not.(ii.eq.i.and.jj.eq.j))then | |
2264 | if(ii.ge.1.and.ii.le.mxgrid | |
2265 | . .and.jj.ge.1.and.jj.le.mxgrid)then | |
2266 | jg=jdropgrid(ii,jj,k) | |
2267 | if(jg.gt.0)then | |
2268 | if(nseg(jg).gt.50)then | |
2269 | x=xptl(n) | |
2270 | y=yptl(n) | |
2271 | vrad=( x*pptl(1,n)/pptl(4,n)+y*pptl(2,n)/pptl(4,n)) | |
2272 | if(vrad.lt.0.)then | |
2273 | ihit=1 | |
2274 | goto94 | |
2275 | endif | |
2276 | endif | |
2277 | endif | |
2278 | endif | |
2279 | endif | |
2280 | enddo | |
2281 | enddo | |
2282 | endif | |
2283 | endif | |
2284 | 94 continue | |
2285 | if(ihit.eq.1)then | |
2286 | istptl(n) = 3 | |
2287 | nseg(jg)=nseg(jg)+1 | |
2288 | delx=2*xgrid/mxgrid*(ii-i) | |
2289 | dely=2*xgrid/mxgrid*(jj-j) | |
2290 | xptl(n)=xptl(n)+delx | |
2291 | yptl(n)=yptl(n)+dely | |
2292 | ! if(k.ge.10.and.k.le.13) | |
2293 | ! . print*,k, ' ',i,j,' ',ii,jj,nseg(jg),' ',delx,dely | |
2294 | endif | |
2295 | 93 continue | |
2296 | endif | |
2297 | ||
2298 | c...prepare /cptl/ for clusters | |
2299 | ||
2300 | if(ish.ge.6)write(ifch,*)'prepare /cptl/ for clusters' | |
2301 | do jj=1,jjj | |
2302 | nptl=nptl+1 | |
2303 | istptl(nptl)=12 | |
2304 | do l=1,4 | |
2305 | pptl(l,nptl)=0. | |
2306 | xorptl(l,nptl)=0 | |
2307 | enddo | |
2308 | sptl(nptl)=0 | |
2309 | uptl(nptl)=0 | |
2310 | optl(nptl)=0 | |
2311 | desptl(nptl)=0 | |
2312 | c limit the maximum number of subcluster to half the number of particle | |
2313 | c (not to have empty subclusters) | |
2314 | nsegmx(jj)=max(1,nint(float(nseg(jj))/float(nsegsuj))) | |
2315 | enddo | |
2316 | ||
2317 | c...prepare /cptl/ for subclusters | |
2318 | ||
2319 | if(ish.ge.6)write(ifch,*)'prepare /cptl/ for subclusters' | |
2320 | mm=0 | |
2321 | do jj=1,jjj | |
2322 | do ii=1,nsegmx(jj) | |
2323 | mm=mm+1 | |
2324 | if(mm.gt.mxcl)stop'jintpo: mxcl too small. ' | |
2325 | mmji(jj,ii)=mm | |
2326 | mseg(mm)=0 | |
2327 | nptl=nptl+1 | |
2328 | istptl(nptl)=10 | |
2329 | do l=1,4 | |
2330 | pptld(l,mm)=0d0 | |
2331 | pptl(l,nptl)=0. | |
2332 | xorptl(l,nptl)=0 | |
2333 | enddo | |
2334 | sptl(nptl)=0 | |
2335 | uptl(nptl)=0 | |
2336 | optl(nptl)=0 | |
2337 | desptl(nptl)=0 | |
2338 | do l=1,nflav | |
2339 | jccl(mm,l,1)=0 | |
2340 | jccl(mm,l,2)=0 | |
2341 | enddo | |
2342 | iorptl(nptl)=nptla+jj | |
2343 | jorptl(n)=0 | |
2344 | if(ii.eq.1)ifrptl(1,nptla+jj)=nptlb+mm | |
2345 | ifrptl(2,nptla+jj)=nptlb+mm | |
2346 | enddo | |
2347 | enddo | |
2348 | ||
2349 | c...separate string segments, add dense area segments to clusters | |
2350 | ||
2351 | if(ish.ge.6)write(ifch,*)'separate string segments' | |
2352 | do 98 n=1,nptla | |
2353 | if(istptl(n).ne.3)goto 98 | |
2354 | i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2355 | j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid) | |
2356 | k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid) | |
2357 | if( i.ge.1.and.i.le.mxgrid | |
2358 | & .and.j.ge.1.and.j.le.mxgrid | |
2359 | & .and.k.ge.1.and.k.le.msgrid)then | |
2360 | jj=jdropgrid(i,j,k) | |
2361 | if(jj.gt.0)then | |
2362 | iimx=nsegmx(jj) | |
2363 | do ii=1,iimx | |
2364 | mm=mmji(jj,ii) | |
2365 | if(mseg(mm).eq.0)goto 10 !not to have an empty cluster | |
2366 | enddo | |
2367 | ii=1+rangen()*iimx | |
2368 | ii=min(ii,iimx) | |
2369 | iini=ii | |
2370 | am2tmpmx=1e30 | |
2371 | 9 ntmp=mmji(jj,ii) | |
2372 | am2tmp=(pptld(4,ntmp)+pptld(3,ntmp)) | |
2373 | & *(pptld(4,ntmp)-pptld(3,ntmp)) | |
2374 | if(am2tmp.gt.5.e3)then | |
2375 | if(am2tmp.lt.am2tmpmx)then | |
2376 | mm=mmji(jj,ii) | |
2377 | am2tmpmx=am2tmp | |
2378 | endif | |
2379 | ii=ii+1 | |
2380 | if(ii.gt.iimx)ii=1 | |
2381 | if(ii.ne.iini)then | |
2382 | goto 9 | |
2383 | else | |
2384 | goto 10 | |
2385 | endif | |
2386 | endif | |
2387 | mm=mmji(jj,ii) | |
2388 | 10 continue | |
2389 | mseg(mm)=mseg(mm)+1 | |
2390 | ifrptl(1,n)=mm !local use of ifrptl | |
2391 | c write(ifch,*)'end',n | |
2392 | c & ,pptl(4,n),pptl(3,n),idptl(n),i,j,k,sptl(n) | |
2393 | p4tmp=0d0 | |
2394 | do l=1,3 | |
2395 | pptld(l,mm)= pptld(l,mm) + dble(pptl(l,n)) | |
2396 | p4tmp=p4tmp+dble(pptl(l,n))*dble(pptl(l,n)) | |
2397 | enddo | |
2398 | p4tmp=sqrt(p4tmp+dble(pptl(5,n)*pptl(5,n))) | |
2399 | pptld(4,mm)= pptld(4,mm) + p4tmp | |
2400 | c if(mm.eq.86)write(ifch,*)'other',n | |
2401 | c & ,pptl(4,n),pptl(3,n),pptl(5,n),idptl(n),k,p4tmp | |
2402 | c & ,pptld(4,mm),pptld(3,mm),pptld(2,mm),pptld(1,mm) | |
2403 | c & ,(pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm)) | |
2404 | xorptl(1,nptlb+mm)=xorptl(1,nptlb+mm)+xptl(n) | |
2405 | xorptl(2,nptlb+mm)=xorptl(2,nptlb+mm)+yptl(n) | |
2406 | xorptl(3,nptlb+mm)=xorptl(3,nptlb+mm)+zptl(n) | |
2407 | xorptl(4,nptlb+mm)=xorptl(4,nptlb+mm)+tptl(n) | |
2408 | sptl(nptlb+mm)=sptl(nptlb+mm)+sptl(n) | |
2409 | aa=cos(phievt) | |
2410 | bb=sin(phievt) | |
2411 | cc=-sin(phievt) | |
2412 | dd=cos(phievt) | |
2413 | xrot=xptl(n)*aa+yptl(n)*bb | |
2414 | yrot=xptl(n)*cc+yptl(n)*dd | |
2415 | uptl(nptlb+mm)=uptl(nptlb+mm)+xrot**2 | |
2416 | optl(nptlb+mm)=optl(nptlb+mm)+yrot**2 | |
2417 | desptl(nptlb+mm)=desptl(nptlb+mm)+xrot*yrot | |
2418 | id=idptl(n) | |
2419 | ida=iabs(id/10) | |
2420 | ids=id/iabs(id) | |
2421 | if(ida.ne.111.and.ida.ne.222.and.ida.ne.333)id=id/10*10 | |
2422 | if(ida.eq.111.or. ida.eq.222.or. ida.eq.333)id=id/10*10+ids | |
2423 | if(ida.eq.213)id=1230*ids | |
2424 | ic(1)=idtrai(1,id,1) | |
2425 | ic(2)=idtrai(2,id,1) | |
2426 | call iddeco(ic,jc) | |
2427 | do l=1,nflav | |
2428 | jccl(mm,l,1)=jccl(mm,l,1)+jc(l,1) | |
2429 | jccl(mm,l,2)=jccl(mm,l,2)+jc(l,2) | |
2430 | enddo | |
2431 | else | |
2432 | idropgrid(i,j,k)=0 | |
2433 | endif | |
2434 | endif | |
2435 | 98 continue | |
2436 | ||
2437 | c...associate segments to clusters | |
2438 | ||
2439 | if(ish.ge.6)write(ifch,*)'associate segments to clusters' | |
2440 | naseg(0)=0 | |
2441 | do jj=1,jjj | |
2442 | do ii=1,nsegmx(jj) | |
2443 | mm=mmji(jj,ii) | |
2444 | naseg(mm)=naseg(mm-1)+mseg(mm) | |
2445 | nfseg(mm)=0 | |
2446 | enddo | |
2447 | enddo | |
2448 | do 97 n=1,nptla | |
2449 | if(istptl(n).ne.3)goto97 | |
2450 | istptl(n) = 2 | |
2451 | mm=ifrptl(1,n) | |
2452 | nfseg(mm)=nfseg(mm)+1 | |
2453 | nsegmt(naseg(mm-1)+nfseg(mm))=n | |
2454 | 97 enddo | |
2455 | do jj=1,jjj | |
2456 | nst=0 | |
2457 | do ii=1,nsegmx(jj) | |
2458 | mm=mmji(jj,ii) | |
2459 | if(mseg(mm).ne.nfseg(mm))stop'jintpo: mseg.ne.nfseg ' | |
2460 | nst=nst+mseg(mm) | |
2461 | enddo | |
2462 | if(nst.ne.nseg(jj))stop'sum(mseg(mm)).ne.nseg(jj)' | |
2463 | enddo | |
2464 | ||
2465 | c...finish cluster storage to /cptl/ | |
2466 | ||
2467 | if(ish.ge.6)write(ifch,*)'finish cluster storage to /cptl/' | |
2468 | xx=0. | |
2469 | yy=0. | |
2470 | xy=0. | |
2471 | mjjsegsum=0 | |
2472 | do jj=1,jjj | |
2473 | njj=nptla+jj | |
2474 | mjjseg=0 | |
2475 | do l=1,nflav | |
2476 | jcjj(l,1)=0 | |
2477 | jcjj(l,2)=0 | |
2478 | enddo | |
2479 | do ii=1,nsegmx(jj) | |
2480 | mm=mmji(jj,ii) | |
2481 | n=nptlb+mm | |
2482 | ||
2483 | do l=1,nflav | |
2484 | jc(l,1)=jccl(mm,l,1) | |
2485 | jc(l,2)=jccl(mm,l,2) | |
2486 | ke(l)=jc(l,1)-jc(l,2) | |
2487 | jcjj(l,1)=jcjj(l,1)+jc(l,1) | |
2488 | jcjj(l,2)=jcjj(l,2)+jc(l,2) | |
2489 | enddo | |
2490 | call idenct(jc,idptl(n) | |
2491 | * ,ibptl(1,n),ibptl(2,n),ibptl(3,n),ibptl(4,n)) | |
2492 | ttest=0d0 | |
2493 | do ji=1,4 | |
2494 | ptest(ji)=0d0 | |
2495 | do ns=naseg(mm-1)+1,naseg(mm) | |
2496 | ni=nsegmt(ns) | |
2497 | ptest(ji)=ptest(ji)+dble(pptl(ji,ni)) | |
2498 | enddo | |
2499 | ptest(ji)=abs(ptest(ji)-pptld(ji,mm)) | |
2500 | ttest=ttest+ptest(ji) | |
2501 | enddo | |
2502 | p52=(pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm)) | |
2503 | & -pptld(1,mm)**2-pptld(2,mm)**2 | |
2504 | pptld(5,mm)=0d0 | |
2505 | if(p52.gt.0d0)then | |
2506 | pptld(5,mm)=sqrt(p52) | |
2507 | jerr(2)=jerr(2)+1 | |
2508 | elseif(p52.le.0d0)then | |
2509 | jerr(3)=jerr(3)+1 | |
2510 | pptld(5,mm)=dble(utamnu(ke(1),ke(2),ke(3),ke(4),ke(5),ke(6),5)) | |
2511 | pptld(4,mm)=sqrt(pptld(3,mm)*pptld(3,mm) | |
2512 | & +pptld(2,mm)*pptld(2,mm) | |
2513 | & +pptld(1,mm)*pptld(1,mm) | |
2514 | & +pptld(5,mm)*pptld(5,mm)) | |
2515 | endif | |
2516 | if(ish.ge.1.and.(abs(ttest).gt.1.d0.or.pptld(5,mm).gt.120d0))then | |
2517 | call utmsg('jintpo&') | |
2518 | write(ifmt,*)'***** Warning in jintpo !',ntry | |
2519 | write(ifch,*)'***** jintpo: momenta messed up (ttest > 0)' | |
2520 | write(ifch,*)'*****',mm,n,mseg(mm),p52,ttest | |
2521 | write(ifch,*)'*****',jj,nsegp4(jj),p4mean(jj) | |
2522 | write(ifch,'(a,10x,4f15.4)')'*****',(pptld(ji,mm),ji=1,4) | |
2523 | do ns=naseg(mm-1)+1,naseg(mm) | |
2524 | ni=nsegmt(ns) | |
2525 | write(ifch,'(a,2i5,5f15.4,f12.4)')'*****',ni,idptl(ni) | |
2526 | * ,(pptl(ji,ni),ji=1,4),pptl(5,ni)**2 | |
2527 | * ,(pptl(4,ni)+pptl(3,ni))*(pptl(4,ni)-pptl(3,ni)) | |
2528 | * -pptl(1,ni)**2-pptl(2,ni)**2 | |
2529 | enddo | |
2530 | endif | |
2531 | if(pptld(5,mm).gt.120d0)then | |
2532 | p4max=0. | |
2533 | nh=0 | |
2534 | do ns=naseg(mm-1)+1,naseg(mm) | |
2535 | ni=nsegmt(ns) | |
2536 | if(pptl(4,ni).ge.p4max)then | |
2537 | nh=ni | |
2538 | p4max=pptl(4,ni) | |
2539 | endif | |
2540 | enddo | |
2541 | if(nh.le.0)then | |
2542 | stop'Cannot be in jintpo ...' | |
2543 | else !put back nh as normal particle | |
2544 | iaaptl(nh)=0 | |
2545 | istptl(nh)=0 | |
2546 | endif | |
2547 | if(ish.ge.1) | |
2548 | & write(ifch,*)'***** Redo cluster without heavy particle :' | |
2549 | & ,nh,ntry | |
2550 | goto 8888 | |
2551 | endif | |
2552 | do l=1,5 | |
2553 | pptl(l,n)=sngl(pptld(l,mm)) | |
2554 | enddo | |
2555 | mjjseg=mjjseg+mseg(mm) | |
2556 | do l=1,4 | |
2557 | pptl(l,njj)=pptl(l,njj)+pptl(l,n) | |
2558 | xorptl(l,njj)=xorptl(l,njj)+xorptl(l,n) | |
2559 | xorptl(l,n)=xorptl(l,n)/float(mseg(mm)) | |
2560 | enddo | |
2561 | sptl(njj)=sptl(njj)+sptl(n) | |
2562 | uptl(njj)=uptl(njj)+uptl(n) | |
2563 | optl(njj)=optl(njj)+optl(n) | |
2564 | desptl(njj)=desptl(njj)+desptl(n) | |
2565 | sptl(n)=sptl(n)/float(mseg(mm)) | |
2566 | uptl(n)=uptl(n)/float(mseg(mm)) | |
2567 | optl(n)=optl(n)/float(mseg(mm)) | |
2568 | desptl(n)=desptl(n)/float(mseg(mm)) | |
2569 | radptl(n)=0 | |
2570 | istptl(n)=10 | |
2571 | ifrptl(1,n)=0 | |
2572 | ifrptl(2,n)=0 | |
2573 | tivptl(1,n)=xorptl(4,n) | |
2574 | tivptl(2,n)=ainfin | |
2575 | ityptl(n)=60 | |
2576 | vocri=pptl(5,n)/epscri(ioclude) | |
2577 | vo=max(vocri,vocell) | |
2578 | radptl(n)=(3.*vo/4./pi)**0.3333 | |
2579 | dezptl(n)=0. | |
2580 | enddo | |
2581 | do l=1,4 | |
2582 | xorptl(l,njj)=xorptl(l,njj)/float(mjjseg) | |
2583 | enddo | |
2584 | mjjsegsum=mjjsegsum+mjjseg | |
2585 | xx=xx+uptl(njj) | |
2586 | yy=yy+optl(njj) | |
2587 | xy=xy+desptl(njj) | |
2588 | sptl(njj)=sptl(njj)/float(mjjseg) | |
2589 | uptl(njj)=uptl(njj)/float(mjjseg) | |
2590 | optl(njj)=optl(njj)/float(mjjseg) | |
2591 | desptl(njj)=desptl(njj)/float(mjjseg) | |
2592 | pjj52=(pptl(4,njj)+pptl(3,njj))*(pptl(4,njj)-pptl(3,njj)) | |
2593 | & -pptl(1,njj)**2-pptl(2,njj)**2 | |
2594 | pptl(5,njj)=0 | |
2595 | if(pjj52.gt.0)pptl(5,njj)=sqrt(pjj52) | |
2596 | ityptl(njj)=60 | |
2597 | call idenct(jc,idptl(njj) | |
2598 | * ,ibptl(1,njj),ibptl(2,njj),ibptl(3,njj),ibptl(4,njj)) | |
2599 | !................. | |
2600 | !inertia tensor: | |
2601 | !<y**2> -<x*y> with <x**2>=uptl <y**2>=optl | |
2602 | !-<x*y> <x**2> <xy>=desptl | |
2603 | !................. | |
2604 | enddo | |
2605 | if(mjjsegsum.gt.0)then | |
2606 | xx=xx/float(mjjsegsum) | |
2607 | yy=yy/float(mjjsegsum) | |
2608 | xy=xy/float(mjjsegsum) | |
2609 | dta=0.5*abs(xx-yy) | |
2610 | ev1=(xx+yy)/2+sqrt(dta**2+xy**2) | |
2611 | ev2=(xx+yy)/2-sqrt(dta**2+xy**2) | |
2612 | if(xy.lt.0..and.dta.ne.0.)then | |
2613 | ranphi=0.5*atan(-xy/dta) | |
2614 | elseif(xy.gt.0..and.dta.ne.0.)then | |
2615 | ranphi=-0.5*atan(xy/dta) | |
2616 | else | |
2617 | ranphi=0 | |
2618 | endif | |
2619 | else | |
2620 | ranphi=0 | |
2621 | endif | |
2622 | ||
2623 | c...print | |
2624 | ||
2625 | if(ish.ge.5)then | |
2626 | write(ifch,*)'print' | |
2627 | do k=1,msgrid | |
2628 | write(ifch,*)'k=',k,' jclu=',jclu(k) | |
2629 | do j=mxgrid,1,-1 | |
2630 | write(ifch,'(10i4,3x,10i4)')(idropgrid(i,j,k),i=1,mxgrid) | |
2631 | & ,(jdropgrid(i,j,k),i=1,mxgrid) | |
2632 | enddo | |
2633 | enddo | |
2634 | write(ifch,'(a,a)') | |
2635 | & ' k jj nseg mm mseg n mass' | |
2636 | & ,' s y z t ' | |
2637 | do jj=1,jjj | |
2638 | do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj)) | |
2639 | mm=mmji(jj,ii) | |
2640 | n=nptlb+mm | |
2641 | sg=pptl(3,n)/abs(pptl(3,n)) | |
2642 | tm=sqrt(pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2) | |
2643 | y=sg*alog((pptl(4,n)+sg*pptl(3,n))/tm) | |
2644 | c if(kclu(jj).eq.44)print *,tm,pptl(4,n),pptl(3,n),iorptl(n) | |
2645 | write(ifch,'(2i5,i6,i8,2i6,5f10.3)') | |
2646 | & kclu(jj),jj,nseg(jj),mm,mseg(mm),n,pptl(5,n) | |
2647 | & ,sptl(n),y,xorptl(3,n),xorptl(4,n) | |
2648 | enddo | |
2649 | enddo | |
2650 | endif | |
2651 | ||
2652 | c...compute ranecc | |
2653 | ||
2654 | weiecc=0 | |
2655 | ranecc=0 | |
2656 | do jj=1,jjj | |
2657 | do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj)) | |
2658 | mm=mmji(jj,ii) | |
2659 | ip=nptlb+mm | |
2660 | ipo=iorptl(ip) | |
2661 | xx=uptl(ipo) ! <x**2> | |
2662 | yy=optl(ipo) ! <y**2> | |
2663 | xy=desptl(ipo) ! <x*y> | |
2664 | dta=0.5*abs(xx-yy) | |
2665 | ev1=(xx+yy)/2+sqrt(dta**2+xy**2) | |
2666 | ev2=(xx+yy)/2-sqrt(dta**2+xy**2) | |
2667 | yy=ev1 | |
2668 | xx=ev2 | |
2669 | ecc=(yy-xx)/(yy+xx) | |
2670 | ranecc=ranecc+ecc*pptl(5,ip) | |
2671 | weiecc=weiecc+pptl(5,ip) | |
2672 | enddo | |
2673 | enddo | |
2674 | if(weiecc.gt.0.)then | |
2675 | ranecc=ranecc/weiecc | |
2676 | endif | |
2677 | ||
2678 | c...decay | |
2679 | ||
2680 | if(ish.ge.6)write(ifch,*)'decay' | |
2681 | if(ifrade.eq.0.or.ispherio.gt.0)goto1000 | |
2682 | if(jdecay.eq.0)goto1000 | |
2683 | ctp060829 npl=0 | |
2684 | do jj=1,jjj | |
2685 | do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj)) | |
2686 | mm=mmji(jj,ii) | |
2687 | np=nptlb+mm | |
2688 | nptlij=nptl | |
2689 | iret=1 | |
2690 | c ishb=ish | |
2691 | c ish=0 | |
2692 | call hnbaaa(np,iret) | |
2693 | c iret=1 | |
2694 | c ish=ishb | |
2695 | if(iret.eq.1)then | |
2696 | istptl(np)=istptl(np)+2 | |
2697 | do ns=naseg(mm-1)+1,naseg(mm) | |
2698 | n=nsegmt(ns) | |
2699 | istptl(n) = 0 | |
2700 | enddo | |
2701 | else | |
2702 | istptl(np)=istptl(np)+1 | |
2703 | ifrptl(1,np)=nptlij+1 | |
2704 | ifrptl(2,np)=nptl | |
2705 | do n=nptlij+1,nptl | |
2706 | iorptl(n)=np | |
2707 | jorptl(n)=0 | |
2708 | istptl(n)=0 | |
2709 | ifrptl(1,n)=0 | |
2710 | ifrptl(2,n)=0 | |
2711 | !xorptl(1,n) already defined in hnbaaa | |
2712 | !xorptl(2,n) already defined in hnbaaa | |
2713 | !xorptl(3,n) already defined in hnbaaa | |
2714 | !xorptl(4,n) already defined in hnbaaa | |
2715 | tivptl(1,n)=xorptl(4,n) | |
2716 | call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm) | |
2717 | r=rangen() | |
2718 | tivptl(2,n)=tivptl(1,n)+taugm*(-alog(r)) | |
2719 | !ityptl(n) already defined in hnbaaa | |
2720 | radptl(n)=0. | |
2721 | dezptl(n)=0. | |
2722 | itsptl(n)=0 | |
2723 | rinptl(n)=kclu(jj)-msgrid/2 | |
2724 | enddo | |
2725 | endif | |
2726 | enddo | |
2727 | enddo | |
2728 | ||
2729 | 1000 continue | |
2730 | call utprix('jintpo',ish,ishini,4) | |
2731 | end | |
2732 | ||
2733 | cc----------------------------------------------------------------------- | |
2734 | c subroutine jrad(i,nq,na,jc,rad) | |
2735 | cc----------------------------------------------------------------------- | |
2736 | cc return hadron radius (data taken from huefner and povh) | |
2737 | cc----------------------------------------------------------------------- | |
2738 | c include 'epos.inc' | |
2739 | c integer jc(nflav,2),kc(nflav) | |
2740 | c | |
2741 | c id=iabs(idptl(i)) | |
2742 | c am=pptl(5,i) | |
2743 | c if(id.lt.10000)then | |
2744 | c k=mod(id,10) | |
2745 | c else | |
2746 | c k=1 | |
2747 | c endif | |
2748 | c do l=1,nflav | |
2749 | c kc(l)=iabs(jc(l,1)-jc(l,2)) | |
2750 | c enddo | |
2751 | c | |
2752 | c if(nq.eq.0)then ! mesons | |
2753 | c if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0.and.kc(4).eq.0)then | |
2754 | c if(k.eq.0)then ! flavor singlet pseudoscalar mesons | |
2755 | c if(am.ge.0.000)then | |
2756 | c rad=0.64 ! pi0 | |
2757 | c if(am.ge.0.500)then | |
2758 | c rad=0.60 ! eta | |
2759 | c if(am.ge.0.900)then | |
2760 | c rad=0.40 ! eta prime | |
2761 | c if(am.ge.2.900)then | |
2762 | c rad=0.17 ! eta charm | |
2763 | c endif | |
2764 | c endif | |
2765 | c endif | |
2766 | c else | |
2767 | c write(ifch,*) | |
2768 | c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am | |
2769 | c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6) | |
2770 | c call utstop('jrad: meson radius not defined&') | |
2771 | c endif | |
2772 | c else ! flavor singlet vector mesons | |
2773 | c if(am.ge.0.000)then | |
2774 | c rad=0.72 ! rho,omega | |
2775 | c if(am.ge.1.000)then | |
2776 | c rad=0.46 ! phi | |
2777 | c if(am.ge.3.000)then | |
2778 | c rad=0.20 ! J/psi | |
2779 | c endif | |
2780 | c endif | |
2781 | c else | |
2782 | c write(ifch,*) | |
2783 | c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am | |
2784 | c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6) | |
2785 | c call utstop('jrad: meson radius not defined&') | |
2786 | c endif | |
2787 | c endif | |
2788 | c elseif(kc(3).eq.0.and.kc(4).eq.0)then ! nonstrange, noncharmed | |
2789 | c if(k.eq.0)then | |
2790 | c rad=0.64 ! pi | |
2791 | c else | |
2792 | c rad=0.72 ! resonances | |
2793 | c endif | |
2794 | c elseif(kc(3).ne.0.and.kc(4).eq.0)then ! strange | |
2795 | c if(k.eq.0)then | |
2796 | c rad=0.59 ! kaons | |
2797 | c else | |
2798 | c rad=0.68 ! kaon resonances | |
2799 | c endif | |
2800 | c else ! charmed | |
2801 | c write(ifch,*)'i:',i,' id:',idptl(i) | |
2802 | c call utstop('jrad: radius of meson not defined&') | |
2803 | c endif | |
2804 | c else !baryons | |
2805 | c if(kc(4).gt.0)then ! charmed | |
2806 | c write(ifch,*) | |
2807 | c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am | |
2808 | c write(ifch,*)'i:',i,' id:',idptl(i) | |
2809 | c call utstop('jrad: radius of charmed baryon not defined&') | |
2810 | c elseif(kc(3).eq.0)then ! nonstrange | |
2811 | c if(k.eq.0)then | |
2812 | c rad=0.82 !nucleons | |
2813 | c else | |
2814 | c rad=1.00 !resonances | |
2815 | c endif | |
2816 | c elseif(kc(3).eq.1)then ! strange | |
2817 | c if(k.eq.0)then | |
2818 | c rad=0.76 !lambda, sigma | |
2819 | c else | |
2820 | c rad=0.93 !resonances | |
2821 | c endif | |
2822 | c elseif(kc(3).eq.2)then ! double strange | |
2823 | c if(k.eq.0)then | |
2824 | c rad=0.71 !cascades | |
2825 | c else | |
2826 | c rad=0.87 !resonances | |
2827 | c endif | |
2828 | c elseif(kc(3).ge.3)then ! triple strange | |
2829 | c rad=0.79 !omega | |
2830 | c else | |
2831 | c write(ifch,*) | |
2832 | c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am | |
2833 | c write(ifch,*) | |
2834 | c * 'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6), | |
2835 | c * ' |q-qbar|:',(kc(l),l=1,6) | |
2836 | c call utstop('jrad: should not happen&') | |
2837 | c endif | |
2838 | cc string fragments with |#q|>3 | |
2839 | c if(na.gt.3)then | |
2840 | c a=(na/3.)**(1./3.) | |
2841 | c if(ish.ge.7)then | |
2842 | c call utmsg('jrad ') | |
2843 | c write(ifch,*) | |
2844 | c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am | |
2845 | c write(ifch,*) | |
2846 | c * 'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6), | |
2847 | c * ' |q-qbar|:',(kc(l),l=1,6) | |
2848 | c write(ifch,*)'nq:',nq,' na:',na,' r:',rad,' ar:',a*rad | |
2849 | c call utmsgf | |
2850 | c endif | |
2851 | c rad=rad*a | |
2852 | c endif | |
2853 | c endif | |
2854 | c | |
2855 | c if(ish.ge.7)then | |
2856 | c write(ifch,*) | |
2857 | c * 'i:',i,' id:',idptl(i),' k:',k,' m:',am,' rad:',rad | |
2858 | c write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6) | |
2859 | c endif | |
2860 | c | |
2861 | c return | |
2862 | c end | |
2863 | c | |
2864 | c----------------------------------------------------------------------- | |
2865 | subroutine jresc | |
2866 | c----------------------------------------------------------------------- | |
2867 | include 'epos.inc' | |
2868 | double precision pa(5),pj(5) | |
2869 | integer ipptl(mxptl) | |
2870 | ||
2871 | call utpri('jresc ',ish,ishini,4) | |
2872 | ||
2873 | iret=0 | |
2874 | nptlpt=maproj+matarg | |
2875 | np=0 | |
2876 | do i=nptlpt+1,nptl | |
2877 | if(istptl(i).eq.0 | |
2878 | * .and.idptl(i).lt.10000.and.pptl(5,i).gt.0.01)then | |
2879 | np=np+1 | |
2880 | ipptl(np)=i | |
2881 | endif | |
2882 | enddo | |
2883 | if(np.lt.2)goto1001 | |
2884 | do ii=1,np | |
2885 | i=ipptl(ii) | |
2886 | if(mod(iabs(idptl(i)),10).lt.2)then | |
2887 | call idmass(idptl(i),ami) | |
2888 | dm=abs(ami-pptl(5,i)) | |
2889 | if(dm.gt.0.001)then | |
2890 | ntry=0 | |
2891 | 1 continue | |
2892 | ntry=ntry+1 | |
2893 | 2 jj=1+int(rangen()*np) | |
2894 | j=ipptl(jj) | |
2895 | if(ish.ge.4)write(ifch,*)i,j,istptl(j) | |
2896 | if(j.eq.i)goto2 | |
2897 | if(mod(iabs(idptl(j)),10).lt.2)then | |
2898 | call idmass(idptl(j),amj) | |
2899 | else | |
2900 | amj=pptl(5,j) | |
2901 | endif | |
2902 | do l=1,5 | |
2903 | pa(l)=dble(pptl(l,i)) | |
2904 | pj(l)=dble(pptl(l,j)) | |
2905 | enddo | |
2906 | if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70) | |
2907 | if(ish.ge.4)write(ifch,11)i,idptl(i),'before:',pa,'want:',ami | |
2908 | if(ish.ge.4)write(ifch,11)j,idptl(j),'before:',pj,'want:',amj | |
2909 | call jrescl(pa,dble(ami),pj,dble(amj),iret) | |
2910 | if(iret.eq.1)then | |
2911 | if(ntry.le.50)then | |
2912 | goto1 | |
2913 | else | |
2914 | goto1001 | |
2915 | endif | |
2916 | endif | |
2917 | if(ish.ge.4)write(ifch,11)i,idptl(i),' after:',pa | |
2918 | if(ish.ge.4)write(ifch,11)j,idptl(j),' after:',pj | |
2919 | if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70) | |
2920 | do l=1,5 | |
2921 | pptl(l,i)=sngl(pa(l)) | |
2922 | pptl(l,j)=sngl(pj(l)) | |
2923 | enddo | |
2924 | endif | |
2925 | endif | |
2926 | enddo | |
2927 | 11 format(i5,1x,i5,1x,a,1x,5(d8.2,1x),a,1x,e8.2) | |
2928 | ||
2929 | 1000 continue | |
2930 | call utprix('jresc ',ish,ishini,4) | |
2931 | return | |
2932 | ||
2933 | 1001 continue | |
2934 | if(ish.ge.1)then | |
2935 | write(ifmt,'(a)')'jresc: could not put on shell' | |
2936 | endif | |
2937 | goto1000 | |
2938 | ||
2939 | end | |
2940 | ||
2941 | c----------------------------------------------------------------------- | |
2942 | subroutine jrescl(p1,am1,p2,am2,iret) | |
2943 | c----------------------------------------------------------------------- | |
2944 | c rescale momenta of two particles such that the masses assume given | |
2945 | c values. | |
2946 | c input: | |
2947 | c p1, p2: momenta of the two particles | |
2948 | c am1, am2: desired masses of the two particles | |
2949 | c output: | |
2950 | c p1, p2: new momenta of the two particles | |
2951 | c----------------------------------------------------------------------- | |
2952 | include 'epos.inc' | |
2953 | double precision p1(5),p2(5) | |
2954 | * ,p1n(5),p2n(5) | |
2955 | * ,a1,a2,a12,am1,am2 | |
2956 | * ,b1,b2,c,d,e,f,g,p,q,r | |
2957 | ||
2958 | call utpri('jrescl',ish,ishini,7) | |
2959 | ||
2960 | iret=0 | |
2961 | a1=p1(5)**2 | |
2962 | a2=p2(5)**2 | |
2963 | a12=p1(4)*p2(4)-p1(3)*p2(3)-p1(2)*p2(2)-p1(1)*p2(1) | |
2964 | if(a12.le.(a1+a2))then | |
2965 | if(ish.ge.7)write(ifch,*)'a_12 < a_1 + a_2' | |
2966 | if(ish.ge.7)write(ifch,*)a12,' < ',a1+a2 | |
2967 | c goto1001 | |
2968 | endif | |
2969 | ||
2970 | 11 format(5(d9.3,1x)) | |
2971 | if(ish.ge.7)write(ifch,11)p1,a1 | |
2972 | if(ish.ge.7)write(ifch,11)p2,a2 | |
2973 | if(ish.ge.7)write(ifch,*)a12 | |
2974 | ||
2975 | c=(a1+a12)/(a2+a12) | |
2976 | d=(a1-am1**2-a2+am2**2)/(a2+a12)*0.5d0 | |
2977 | ||
2978 | e=a1-2d0*a12*c+a2*c**2 | |
2979 | f=2d0*(a1-a12*(c+d)+a2*c*d) | |
2980 | g=a1-2d0*a12*d+a2*d**2-am1**2 | |
2981 | ||
2982 | p=f/e | |
2983 | q=g/e | |
2984 | r=p**2-4d0*q | |
2985 | ||
2986 | if(ish.ge.7)write(ifch,*)'c:',c,' d:',d | |
2987 | if(ish.ge.7)write(ifch,*)'e:',e,' f:',f,' g:',g | |
2988 | if(ish.ge.7)write(ifch,*)'p:',p,' q:',q,' r:',r | |
2989 | if(r.lt.0d0)goto1001 | |
2990 | ||
2991 | b1=-0.5d0*(p-dsqrt(r)) | |
2992 | ||
2993 | b2=b1*c+d | |
2994 | ||
2995 | if(ish.ge.7)write(ifch,*)'b_1:',b1,' b_2:',b2 | |
2996 | ||
2997 | do i=1,4 | |
2998 | p1n(i)=(1d0+b1)*p1(i)-b2*p2(i) | |
2999 | p2n(i)=(1d0+b2)*p2(i)-b1*p1(i) | |
3000 | enddo | |
3001 | ||
3002 | a1=p1n(4)**2-p1n(3)**2-p1n(2)**2-p1n(1)**2 | |
3003 | a2=p2n(4)**2-p2n(3)**2-p2n(2)**2-p2n(1)**2 | |
3004 | if(a1.gt.0d0.and.a2.gt.0d0)then | |
3005 | do i=1,4 | |
3006 | p1(i)=p1n(i) | |
3007 | p2(i)=p2n(i) | |
3008 | enddo | |
3009 | p1(5)=dsqrt(a1) | |
3010 | p2(5)=dsqrt(a2) | |
3011 | if(ish.ge.7)write(ifch,11)p1,a1 | |
3012 | if(ish.ge.7)write(ifch,11)p2,a2 | |
3013 | else | |
3014 | goto1001 | |
3015 | endif | |
3016 | ||
3017 | if(p1(4).lt.0..or.p2(4).lt.0.)goto1001 | |
3018 | ||
3019 | 1000 continue | |
3020 | call utprix('jrescl',ish,ishini,7) | |
3021 | return | |
3022 | ||
3023 | 1001 continue | |
3024 | iret=1 | |
3025 | goto1000 | |
3026 | end | |
3027 | ||
3028 | c----------------------------------------------------------------------- | |
3029 | subroutine jtain(i,x,y,z,t,n,iopt) | |
3030 | c----------------------------------------------------------------------- | |
3031 | c returns intersection (x,y,z,t) of ptl-i-trajectory with taus-line. | |
3032 | c input: | |
3033 | c i: particle number | |
3034 | c iopt: formation time considered (0) or not (1) | |
3035 | c output: | |
3036 | c x,y,z,t: 4-vector of intersection point | |
3037 | c n: exit code | |
3038 | c n=0: ok | |
3039 | c n=1: ptl lives later | |
3040 | c n=2: ptl lives earlier | |
3041 | c n=9: tiv1>tiv2 | |
3042 | c----------------------------------------------------------------------- | |
3043 | include 'epos.inc' | |
3044 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3045 | common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3046 | double precision vv,zza,zz,tt,xo3,xo4,ti1,ti2,derr,dd | |
3047 | double precision ttp,zzp,ttt,zzt,vvt,vvp,spt2m2E,p4 | |
3048 | common/ctfi/tin,tfi | |
3049 | double precision ttau0 | |
3050 | common/cttau0/ttau0 | |
3051 | ||
3052 | n=0 | |
3053 | ||
3054 | tin=0 | |
3055 | tfi=0 | |
3056 | ||
3057 | derr=1d-2 | |
3058 | ttp=tpro*ttaus | |
3059 | zzp=zpro*ttaus | |
3060 | ttt=ttar*ttaus | |
3061 | zzt=ztar*ttaus | |
3062 | vv=sign(min(1.d0,abs(dble(pptl(3,i)))/dble(pptl(4,i))) | |
3063 | & ,dble(pptl(3,i))) | |
3064 | ||
3065 | ||
3066 | if(abs(vv).ge.1.d0)then | |
3067 | spt2m2E=dble(pptl(1,i)*pptl(1,i)+pptl(2,i)*pptl(2,i) | |
3068 | & +pptl(5,i)*pptl(5,i)) | |
3069 | c if(pptl(4,i).le.0.)then | |
3070 | p4=sqrt(dble(pptl(3,i)*pptl(3,i))+spt2m2E) | |
3071 | c else | |
3072 | c p4=dble(pptl(4,i)) | |
3073 | c endif | |
3074 | ctp to avoid precision problem, replace abs(p3)/p4 by sqrt(1-(pt2+m2)/E2) | |
3075 | spt2m2E=min(1.d0,sqrt(spt2m2E)/p4) | |
3076 | vv=sign(sqrt((1d0+spt2m2E)*(1d0-spt2m2E)),dble(pptl(3,i))) | |
3077 | endif | |
3078 | xo3=dble(xorptl(3,i)) | |
3079 | xo4=dble(xorptl(4,i)) | |
3080 | zza=xo3-xo4*vv | |
3081 | if(iopt.eq.0)ti1=dble(tivptl(1,i)) | |
3082 | if(iopt.eq.1)ti1=dble(xo4) | |
3083 | ti2=dble(tivptl(2,i)) | |
3084 | ||
3085 | if(ti1.gt.ti2)then | |
3086 | n=9 | |
3087 | goto1 | |
3088 | endif | |
3089 | ||
3090 | zfi=sngl(xo3+(ti2-xo4)*vv) | |
3091 | call jtaus(zfi,tzfi,szfi) | |
3092 | tfi=tzfi | |
3093 | if(tfi.ge.sngl(ti2))then | |
3094 | n=2 | |
3095 | goto1 | |
3096 | endif | |
3097 | ||
3098 | zin=sngl(xo3+(ti1-xo4)*vv) | |
3099 | call jtaus(zin,tzin,szin) | |
3100 | tin=tzin | |
3101 | if(tin.le.sngl(ti1))then | |
3102 | n=1 | |
3103 | goto1 | |
3104 | endif | |
3105 | ||
3106 | ||
3107 | 1 continue | |
3108 | ||
3109 | if(ttaus.le.ttau0)then | |
3110 | tt=ttaus | |
3111 | zz=xo3+(tt-xo4)*vv | |
3112 | if(tt.lt.ti1.and.n.eq.0)n=1 | |
3113 | if(tt.ge.ti2.and.n.eq.0)n=2 | |
3114 | goto1000 | |
3115 | else | |
3116 | vvt=zzt/ttt | |
3117 | vvp=zzp/ttp | |
3118 | tt=(ttt+(zza-zzt)*vvt)/(1-vv*vvt) | |
3119 | zz=xo3+(tt-xo4)*vv | |
3120 | if(zz.le.zzt)then | |
3121 | if(tt.lt.ti1.and.n.eq.0)n=1 | |
3122 | if(tt.ge.ti2.and.n.eq.0)n=2 | |
3123 | goto1000 | |
3124 | endif | |
3125 | tt=(ttp+(zza-zzp)*vvp)/(1-vv*vvp) | |
3126 | zz=xo3+(tt-xo4)*vv | |
3127 | if(zz.ge.zzp)then | |
3128 | if(tt.lt.ti1.and.n.eq.0)n=1 | |
3129 | if(tt.ge.ti2.and.n.eq.0)n=2 | |
3130 | goto1000 | |
3131 | endif | |
3132 | dd=1-vv**2 | |
3133 | if(sngl(dd).eq.0..and.vv.gt.0.)then | |
3134 | tt=-(ttaus**2+zza**2)/2d0/zza | |
3135 | elseif(sngl(dd).eq.0..and.vv.lt.0.)then | |
3136 | tt=(ttaus**2+zza**2)/2d0/zza | |
3137 | else | |
3138 | tt=(zza*vv+dsqrt(zza**2+ttaus**2*dd))/dd | |
3139 | endif | |
3140 | zz=xo3+(tt-xo4)*vv | |
3141 | if(tt.lt.ti1.and.n.eq.0)n=1 | |
3142 | if(tt.ge.ti2.and.n.eq.0)n=2 | |
3143 | if(dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr*ttaus**2.and. | |
3144 | *dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr)then | |
3145 | if(ish.ge.1)then | |
3146 | call utmsg('jtain') | |
3147 | write(ifch,*)'***** ttaus**2 .ne. (tt+zz)*(tt-zz)' | |
3148 | write(ifch,*)sngl(ttaus**2),sngl((tt+zz)*(tt-zz)) | |
3149 | call utmsgf | |
3150 | endif | |
3151 | goto1000 | |
3152 | endif | |
3153 | endif | |
3154 | ||
3155 | 1000 t=sngl(tt) | |
3156 | z=sngl(zz) | |
3157 | x=xorptl(1,i)+(t-xorptl(4,i))*pptl(1,i)/pptl(4,i) | |
3158 | y=xorptl(2,i)+(t-xorptl(4,i))*pptl(2,i)/pptl(4,i) | |
3159 | return | |
3160 | end | |
3161 | ||
3162 | c----------------------------------------------------------------------- | |
3163 | subroutine jtaix(i,tau,zor,tor,z,t) | |
3164 | c----------------------------------------------------------------------- | |
3165 | c returns intersection z,t of ptl-i-trajectory with hyperbola h. | |
3166 | c h: (t-tor)**2-(z-zor)**2=tau**2 . | |
3167 | c zor, tor double precision. | |
3168 | c----------------------------------------------------------------------- | |
3169 | include 'epos.inc' | |
3170 | double precision tor,zor,tors,zors,vv,cc,dd,ttau,derr,tt,zz | |
3171 | derr=1d-3 | |
3172 | ttau=dble(tau) | |
3173 | zors=dble(xorptl(3,i))-zor | |
3174 | tors=dble(xorptl(4,i))-tor | |
3175 | vv=dble(pptl(3,i)/pptl(4,i)) | |
3176 | vv=dmin1(vv,1d0) | |
3177 | vv=dmax1(vv,-1d0) | |
3178 | cc=zors-tors*vv | |
3179 | dd=1d0-vv**2 | |
3180 | dd=dmax1(dd,0d0) | |
3181 | if(dd.eq.0d0.and.cc.eq.0d0)then | |
3182 | if(tau.eq.0.)tt=0d0 | |
3183 | if(tau.ne.0.)tt=dble(ainfin) | |
3184 | zz=tt | |
3185 | goto1000 | |
3186 | elseif(dd.eq.0d0)then | |
3187 | tt=-(ttau**2+cc**2)/2d0/cc/vv | |
3188 | elseif(dd.lt.1e-8)then | |
3189 | tt=-(ttau**2+cc**2)/2d0/cc/vv | |
3190 | call utmsg('jtaix') | |
3191 | write(ifch,*)'***** dd = ',dd,' treated as zero' | |
3192 | call utmsgf | |
3193 | else | |
3194 | tt=(cc*vv+dsqrt(cc**2+ttau**2*dd)) | |
3195 | tt=tt/dd | |
3196 | endif | |
3197 | zz=cc+tt*vv | |
3198 | if(dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr*ttau**2 | |
3199 | *.and.dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr | |
3200 | *.and.tors**2-zors**2.lt.1e6)then | |
3201 | if(ish.ge.1)then | |
3202 | call utmsg('jtaix') | |
3203 | write(ifch,*)'***** ttau**2 .ne. (tt+zz)*(tt-zz)' | |
3204 | write(ifch,*)sngl(ttau**2),sngl((tt+zz)*(tt-zz)) | |
3205 | write(ifch,*)'tau,t,z:' | |
3206 | write(ifch,*)tau,tt,zz | |
3207 | write(ifch,*)'#,id(ptl):',i,idptl(i) | |
3208 | write(ifch,*)'zor,tor(str):',zor,tor | |
3209 | write(ifch,*)'zors,tors,p,e(ptl):' | |
3210 | write(ifch,*)sngl(zors),sngl(tors),pptl(3,i),pptl(4,i) | |
3211 | call utmsgf | |
3212 | endif | |
3213 | endif | |
3214 | 1000 z=sngl(zz+zor) | |
3215 | t=sngl(tt+tor) | |
3216 | return | |
3217 | end | |
3218 | ||
3219 | c----------------------------------------------------------------------- | |
3220 | subroutine jtaug(su,so,g,y) | |
3221 | c----------------------------------------------------------------------- | |
3222 | c returns g factor and rapidity y for given su, so | |
3223 | c----------------------------------------------------------------------- | |
3224 | include 'epos.inc' | |
3225 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3226 | common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3227 | double precision ttp,zzp,ttt,zzt,ssp,sst,ssu,sso,ss1,ss2,gg | |
3228 | *,ssav,yyav,hh | |
3229 | double precision ttau0 | |
3230 | common/cttau0/ttau0 | |
3231 | ||
3232 | ssu=dble(su) | |
3233 | sso=dble(so) | |
3234 | ||
3235 | if(ssu.ge.sso)then | |
3236 | sso=(ssu+sso)*0.5d0 + dble(abs(dezzer))*ttaus*0.5d0 | |
3237 | ssu=(ssu+sso)*0.5d0 - dble(abs(dezzer))*ttaus*0.5d0 | |
3238 | so=real(sso) | |
3239 | su=real(ssu) | |
3240 | endif | |
3241 | if(ssu.ge.sso)then | |
3242 | print*,ssu,sso,dble(abs(dezzer))*ttaus*0.5d0 | |
3243 | stop'STOP: sr jtaug: ssu.ge.sso' | |
3244 | endif | |
3245 | ||
3246 | g=1 | |
3247 | ||
3248 | if(ttaus.le.ttau0)return | |
3249 | ||
3250 | ttp=tpro*ttaus | |
3251 | zzp=zpro*ttaus | |
3252 | ttt=ttar*ttaus | |
3253 | zzt=ztar*ttaus | |
3254 | ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp)) | |
3255 | sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt)) | |
3256 | ||
3257 | ssav=(ssu+sso)/2d0 | |
3258 | yyav=ssav/ttaus | |
3259 | if(ssav.ge.ssp)yyav=detap | |
3260 | if(ssav.le.sst)yyav=detat | |
3261 | ||
3262 | gg=0 | |
3263 | if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu) | |
3264 | if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu)) | |
3265 | if(ssu.lt.ssp.and.sso.gt.sst)then | |
3266 | ss1=dmax1(ssu,sst) | |
3267 | ss2=dmin1(sso,ssp) | |
3268 | gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) ) | |
3269 | endif | |
3270 | gg=gg/(sso-ssu) | |
3271 | ||
3272 | hh=0 | |
3273 | if(ssu.lt.sst)hh=hh + dsinh(detat-yyav) * (dmin1(sst,sso)-ssu) | |
3274 | if(sso.gt.ssp)hh=hh + dsinh(detap-yyav) * (sso-dmax1(ssp,ssu)) | |
3275 | if(ssu.lt.ssp.and.sso.gt.sst)then | |
3276 | ss1=dmax1(ssu,sst) | |
3277 | ss2=dmin1(sso,ssp) | |
3278 | hh=hh+ttaus*( dcosh(ss2/ttaus-yyav)-dcosh(ss1/ttaus-yyav) ) | |
3279 | endif | |
3280 | hh=hh/(sso-ssu) | |
3281 | ||
3282 | yyav=yyav+0.5d0*dlog((gg+hh)/(gg-hh)) | |
3283 | ||
3284 | gg=0 | |
3285 | if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu) | |
3286 | if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu)) | |
3287 | if(ssu.lt.ssp.and.sso.gt.sst)then | |
3288 | ss1=dmax1(ssu,sst) | |
3289 | ss2=dmin1(sso,ssp) | |
3290 | gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) ) | |
3291 | endif | |
3292 | gg=gg/(sso-ssu) | |
3293 | ||
3294 | g=sngl(gg) | |
3295 | y=sngl(yyav) | |
3296 | return | |
3297 | end | |
3298 | ||
3299 | c----------------------------------------------------------------------- | |
3300 | subroutine jtaui(s,ts,zs) | |
3301 | c----------------------------------------------------------------------- | |
3302 | c returns time ts and coord zs corresponding to ttaus and inv. length s | |
3303 | c----------------------------------------------------------------------- | |
3304 | ||
3305 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3306 | common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3307 | double precision ttau0 | |
3308 | common/cttau0/ttau0 | |
3309 | ||
3310 | double precision ttp,zzp,ttt,zzt,ssp,sst,ss,zeta | |
3311 | ||
3312 | zs=s | |
3313 | ts=sngl(ttaus) | |
3314 | ||
3315 | if(ttaus.le.ttau0)return | |
3316 | ||
3317 | ttp=tpro*ttaus | |
3318 | zzp=zpro*ttaus | |
3319 | ttt=ttar*ttaus | |
3320 | zzt=ztar*ttaus | |
3321 | ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp)) | |
3322 | sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt)) | |
3323 | ss=dble(s) | |
3324 | ||
3325 | if(ss.le.sst)then | |
3326 | zs=sngl(zzt+ttar*(ss-sst)) | |
3327 | ts=sngl(ttt+(dble(zs)-zzt)*zzt/ttt) | |
3328 | elseif(ss.ge.ssp)then | |
3329 | zs=sngl(zzp+tpro*(ss-ssp)) | |
3330 | ts=sngl(ttp+(dble(zs)-zzp)*zzp/ttp) | |
3331 | else | |
3332 | zeta=ss/ttaus | |
3333 | ts=sngl(ttaus*dcosh(zeta)) | |
3334 | zs=sngl(ttaus*dsinh(zeta)) | |
3335 | endif | |
3336 | ||
3337 | return | |
3338 | end | |
3339 | ||
3340 | c----------------------------------------------------------------------- | |
3341 | subroutine jtauin | |
3342 | c----------------------------------------------------------------------- | |
3343 | c initializes equal time hyperbola at ttaus | |
3344 | c----------------------------------------------------------------------- | |
3345 | include 'epos.inc' | |
3346 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3347 | common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3348 | double precision ttau0,rcproj,rctarg | |
3349 | common/geom1/rcproj,rctarg | |
3350 | common/cttau0/ttau0 | |
3351 | ||
3352 | call utpri('jtauin',ish,ishini,6) | |
3353 | ||
3354 | if(ttaus.gt.ttau0)then | |
3355 | if(rcproj.gt.1d-10)then | |
3356 | detap=dmin1(dble((ypjtl-yhaha)*etafac),dlog(ttaus/rcproj)) | |
3357 | else | |
3358 | detap=dble((ypjtl-yhaha)*etafac) | |
3359 | endif | |
3360 | if(rctarg.gt.1d-10)then | |
3361 | detat=dmax1(dble(-yhaha*etafac),dlog(rctarg/ttaus)) | |
3362 | else | |
3363 | detat=dble(-yhaha*etafac) | |
3364 | endif | |
3365 | tpro=dcosh(detap) | |
3366 | zpro=dsinh(detap) | |
3367 | ttar=dcosh(detat) | |
3368 | ztar=dsinh(detat) | |
3369 | else | |
3370 | detap=0d0 | |
3371 | detat=0d0 | |
3372 | tpro=0d0 | |
3373 | zpro=0d0 | |
3374 | ttar=0d0 | |
3375 | ztar=0d0 | |
3376 | endif | |
3377 | ||
3378 | if(ish.ge.6)then | |
3379 | write(ifch,*)'hyperbola at tau=',ttaus | |
3380 | write(ifch,*)'r_p:',rcproj,' r_t:',rctarg | |
3381 | write(ifch,*)'y_p:',detap,' y_t:',detat | |
3382 | write(ifch,*)'t_p:',tpro,' z_p:',zpro | |
3383 | write(ifch,*)'t_t:',ttar,' z_t:',ztar | |
3384 | endif | |
3385 | ||
3386 | call utprix('jtauin',ish,ishini,6) | |
3387 | return | |
3388 | end | |
3389 | ||
3390 | c----------------------------------------------------------------------- | |
3391 | subroutine jtaus(z,tz,sz) | |
3392 | c----------------------------------------------------------------------- | |
3393 | c returns time tz and inv length sz corresponding to ttaus and z | |
3394 | c----------------------------------------------------------------------- | |
3395 | include 'epos.inc' | |
3396 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3397 | common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3398 | double precision ttau0 | |
3399 | common/cttau0/ttau0 | |
3400 | ||
3401 | double precision ttp,zzp,ttt,zzt,zz,tzz | |
3402 | ||
3403 | tz=sngl(ttaus) | |
3404 | sz=z | |
3405 | ||
3406 | if(ttaus.le.ttau0)return | |
3407 | ||
3408 | ttp=tpro*ttaus | |
3409 | zzp=zpro*ttaus | |
3410 | ttt=ttar*ttaus | |
3411 | zzt=ztar*ttaus | |
3412 | zz=dble(z) | |
3413 | ||
3414 | if(zz.le.zzt)then | |
3415 | tz=sngl(ttt+(zz-zzt)*zzt/ttt) | |
3416 | sz=sngl(ttaus*detat+(zz-zzt)/ttar) | |
3417 | elseif(zz.ge.zzp)then | |
3418 | tz=sngl(ttp+(zz-zzp)*zzp/ttp) | |
3419 | sz=sngl(ttaus*detap+(zz-zzp)/tpro) | |
3420 | else | |
3421 | if(sngl(ttaus).ge.ainfin)then | |
3422 | tz=sngl(ttaus) | |
3423 | sz=0. | |
3424 | if(ish.ge.1)then | |
3425 | call utmsg('jtaus') | |
3426 | write(ifch,*)'***** large ttaus; set tz=ttaus, sz=0' | |
3427 | write(ifch,*)'ttaus=',ttaus,'zz=',zz | |
3428 | call utmsgf | |
3429 | endif | |
3430 | else | |
3431 | tzz=dsqrt(ttaus**2+zz**2) | |
3432 | tz=sngl(tzz) | |
3433 | sz=sngl(ttaus*0.5d0*dlog((tzz+zz)/(tzz-zz))) | |
3434 | endif | |
3435 | endif | |
3436 | ||
3437 | return | |
3438 | end | |
3439 | ||
3440 | c----------------------------------------------------------------------- | |
3441 | subroutine jtaux(t,z,ttaux) | |
3442 | c----------------------------------------------------------------------- | |
3443 | c returns ttaux (-> tau-line) for given t and z. | |
3444 | c ttaux: double precision. | |
3445 | c----------------------------------------------------------------------- | |
3446 | double precision ttaux,tt,zz,rcproj,rctarg,zt1,zp1,zt2,zp2,ttau0 | |
3447 | common/geom1/rcproj,rctarg | |
3448 | common/cttau0/ttau0 | |
3449 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3450 | common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3451 | ||
3452 | tt=dble(t) | |
3453 | zz=dble(z) | |
3454 | ||
3455 | if(tt.gt.ttau0)then | |
3456 | zt1=rctarg-tt | |
3457 | zp1=tt-rcproj | |
3458 | zt2=ztar/ttar*tt | |
3459 | zp2=zpro/tpro*tt | |
3460 | if(zz.le.dmax1(zt1,zt2))then | |
3461 | if(zt1.gt.zt2)then | |
3462 | ttaux=rctarg*dsqrt((tt-zz)/(2d0*rctarg-tt-zz)) | |
3463 | else | |
3464 | ttaux=(ttar*tt-ztar*zz)/(ttar**2-ztar**2) | |
3465 | endif | |
3466 | elseif(zz.ge.dmin1(zp1,zp2))then | |
3467 | if(zp1.lt.zp2)then | |
3468 | ttaux=rcproj*dsqrt((tt+zz)/(2d0*rcproj-tt+zz)) | |
3469 | else | |
3470 | ttaux=(tpro*tt-zpro*zz)/(tpro**2-zpro**2) | |
3471 | endif | |
3472 | else | |
3473 | ttaux=dsqrt(tt**2-zz**2) | |
3474 | endif | |
3475 | else | |
3476 | ttaux=tt | |
3477 | endif | |
3478 | ||
3479 | return | |
3480 | end | |
3481 | ||
3482 | c----------------------------------------------------------------------- | |
3483 | subroutine xjden1(ii,itau,x,y,rad,o,u) | |
3484 | c----------------------------------------------------------------------- | |
3485 | c ii=0: initialization | |
3486 | c ii=1: determining dense regions in space of individual events | |
3487 | c x,y,rad: tranverse coordinates and radius of particle i | |
3488 | c o,u: specifies long range: u < s < o (s: long coordinate) | |
3489 | c ii=2: plot of individual event | |
3490 | c xpar1: itau ; valid: 1,..,10 | |
3491 | c xpar2: z-range: -xpar2 < z < xpar2 | |
3492 | c xpar3, x-range: -xpar3 < x < xpar3 | |
3493 | c xpar4, y-range: -xpar4 < y < xpar4 | |
3494 | c----------------------------------------------------------------------- | |
3495 | include "epos.inc" | |
3496 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3497 | common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3498 | ||
3499 | if(idensi.ne.1)stop'STOP in xjden1: idensi must be set 1' | |
3500 | ||
3501 | dlcoox=0.5 | |
3502 | dlcooy=0.5 | |
3503 | ||
3504 | if(ii.eq.0)then | |
3505 | ||
3506 | do i=1,nzeta | |
3507 | do j=1,mxcoox | |
3508 | do k=1,mxcooy | |
3509 | kdensi(itau,i,j,k)=0 | |
3510 | enddo | |
3511 | enddo | |
3512 | enddo | |
3513 | ||
3514 | elseif(ii.eq.1)then | |
3515 | ||
3516 | if(itau.lt.1.or.itau.gt.mxtau)return | |
3517 | ||
3518 | tau=sngl(ttaus) | |
3519 | zu=u/tau | |
3520 | zo=o/tau | |
3521 | ||
3522 | do 1 i=1,nzeta | |
3523 | zi=-nzeta*dlzeta/2-dlzeta/2+i*dlzeta | |
3524 | if(zu.gt.zi.or.zo.lt.zi)goto1 | |
3525 | do 2 j=1,mxcoox | |
3526 | xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox | |
3527 | do 3 k=1,mxcooy | |
3528 | yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy | |
3529 | if((x-xj)**2+(y-yk)**2.gt.rad**2)goto3 | |
3530 | kdensi(itau,i,j,k)=1 | |
3531 | 3 continue | |
3532 | 2 continue | |
3533 | 1 continue | |
3534 | ||
3535 | elseif(ii.eq.2)then | |
3536 | ||
3537 | itaux=nint(xpar1) | |
3538 | if(itaux.gt.mxtau)stop'STOP in xjden1: itaux too large' | |
3539 | ||
3540 | iz=nint(xpar2/dlzeta) | |
3541 | ix=nint(xpar3/dlcoox) | |
3542 | iy=nint(xpar4/dlcooy) | |
3543 | if(iz.gt.nzeta/2)stop'STOP in xjden1: zeta-range too large' | |
3544 | if(ix.gt.mxcoox/2)stop'STOP in xjden1: x-range too large' | |
3545 | if(iy.gt.mxcooy/2)stop'STOP in xjden1: y-range too large' | |
3546 | ||
3547 | do k=mxcooy/2+1-iy,mxcooy/2+iy | |
3548 | write(ifhi,'(a,f7.2)') '! tau: ',tauv(itaux) | |
3549 | write(ifhi,'(a)') 'openhisto' | |
3550 | write(ifhi,'(a,2f7.2)')'xrange',-iz*dlzeta,iz*dlzeta | |
3551 | write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox | |
3552 | write(ifhi,'(a)') 'set ityp2d 3' | |
3553 | write(ifhi,'(a)') 'txt "xaxis space-time rapidity [z]"' | |
3554 | write(ifhi,'(a)') 'txt "yaxis transverse coordinate x (fm)"' | |
3555 | write(ifhi,'(a,i4)') 'array2d',2*iz | |
3556 | do j=mxcoox/2+1-ix,mxcoox/2+ix | |
3557 | write(ifhi,'(40i2)') (kdensi(itaux,i,j,k), | |
3558 | * i=nzeta/2+1-iz,nzeta/2+iz) | |
3559 | enddo | |
3560 | write(ifhi,'(a)') ' endarray' | |
3561 | write(ifhi,'(a)') 'closehisto plot2d' | |
3562 | enddo | |
3563 | ||
3564 | else | |
3565 | ||
3566 | stop'STOP in xjden1: wrong option' | |
3567 | ||
3568 | endif | |
3569 | ||
3570 | return | |
3571 | end | |
3572 | ||
3573 | c----------------------------------------------------------------------- | |
3574 | subroutine xjden2(ii,itau,x,y,rad,s) | |
3575 | c----------------------------------------------------------------------- | |
3576 | c ii=0: initialization | |
3577 | c ii=1: determining dense regions in space of individual events | |
3578 | c x,y,rad: tranverse coordinates and radius of particle i | |
3579 | c s: long coordinate | |
3580 | c ii=2: plot of individual event | |
3581 | c xpar1: itau ; valid: 1,..,10 | |
3582 | c xpar2: s-range: -xpar2 < s < xpar2 | |
3583 | c xpar3, x-range: -xpar3 < x < xpar3 | |
3584 | c xpar4, y-range: -xpar4 < y < xpar4 | |
3585 | c----------------------------------------------------------------------- | |
3586 | include "epos.inc" | |
3587 | double precision tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3588 | common/cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat | |
3589 | parameter (mxcoos=60) | |
3590 | common/cdensh/kdensh(matau,mxcoos,mxcoox,mxcooy),ktot(matau) | |
3591 | character cy*3 | |
3592 | ||
3593 | dlcoox=0.5 | |
3594 | dlcooy=0.5 | |
3595 | dlcoos=0.5 | |
3596 | ||
3597 | if(ii.eq.0)then | |
3598 | ||
3599 | do i=1,mxcoos | |
3600 | do j=1,mxcoox | |
3601 | do k=1,mxcooy | |
3602 | kdensh(itau,i,j,k)=0 | |
3603 | enddo | |
3604 | enddo | |
3605 | enddo | |
3606 | ktot(itau)=0 | |
3607 | ||
3608 | elseif(ii.eq.1)then | |
3609 | ||
3610 | if(itau.lt.1.or.itau.gt.mxtau)return | |
3611 | ||
3612 | tau=sngl(ttaus) | |
3613 | z=s/tau | |
3614 | ||
3615 | do 1 i=1,mxcoos | |
3616 | si=-mxcoos*dlcoos/2-dlcoos/2+i*dlcoos | |
3617 | do 2 j=1,mxcoox | |
3618 | xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox | |
3619 | do 3 k=1,mxcooy | |
3620 | yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy | |
3621 | if(((x-xj)**2+(y-yk)**2+(z-si)**2).gt.rad**2)goto3 | |
3622 | kdensh(itau,i,j,k)=kdensh(itau,i,j,k)+1 | |
3623 | ktot(itau)=ktot(itau)+1 | |
3624 | 3 continue | |
3625 | 2 continue | |
3626 | 1 continue | |
3627 | ||
3628 | elseif(ii.eq.2)then | |
3629 | ||
3630 | itaux=nint(xpar1) | |
3631 | if(itaux.gt.mxtau)stop'STOP in xjden2: itaux too large' | |
3632 | ||
3633 | is=nint(xpar2/dlcoos) | |
3634 | ix=nint(xpar3/dlcoox) | |
3635 | iy=nint(xpar4/dlcooy) | |
3636 | if(is.gt.mxcoos/2)stop'STOP in xjden2: s-range too large' | |
3637 | if(ix.gt.mxcoox/2)stop'STOP in xjden2: x-range too large' | |
3638 | if(iy.gt.mxcooy/2)stop'STOP in xjden2: y-range too large' | |
3639 | ||
3640 | do k=mxcooy/2+1-iy,mxcooy/2+iy | |
3641 | write(cy,'(f3.1)')-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy | |
3642 | write(ifhi,'(a)') 'openhisto' | |
3643 | write(ifhi,'(a,2f7.2)')'xrange',-is*dlcoos,is*dlcoos | |
3644 | write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox | |
3645 | write(ifhi,'(a)') 'set ityp2d 5' | |
3646 | write(ifhi,'(a)') 'txt "xaxis [z] "' | |
3647 | write(ifhi,'(a)') | |
3648 | *'txt "yaxis x (fm), y='//cy//' fm"' | |
3649 | write(ifhi,'(a,i4)') 'array2d',2*is | |
3650 | do j=mxcoox/2+1-ix,mxcoox/2+ix | |
3651 | do i=mxcoos/2+1-is,mxcoos/2+is | |
3652 | write(ifhi,'(e11.3)') | |
3653 | *kdensh(itaux,i,j,k)/dlcooy/dlcoos/dlcoox/ktot(itaux) | |
3654 | enddo | |
3655 | enddo | |
3656 | write(ifhi,'(a)') ' endarray' | |
3657 | write(ifhi,'(a)') 'closehisto plot2d' | |
3658 | enddo | |
3659 | ||
3660 | else | |
3661 | ||
3662 | stop'STOP in xjden2: wrong option' | |
3663 | ||
3664 | endif | |
3665 | ||
3666 | return | |
3667 | end | |
3668 | ||
3669 | cc----------------------------------------------------------------------- | |
3670 | c subroutine postscript(iii,ii,ic) | |
3671 | cc----------------------------------------------------------------------- | |
3672 | c include 'epos.inc' | |
3673 | c character*10 color(10) | |
3674 | c if(iii.eq.0)then | |
3675 | c ifps=21 | |
3676 | c open(unit=ifps,file='zzz.ps',status='unknown') | |
3677 | c WRITE(ifps,'(a)') '%!PS-Adobe-2.0' | |
3678 | c WRITE(ifps,'(a)') '%%Title: tt2.fig' | |
3679 | c WRITE(ifps,'(a)') '%%Orientation: Portrait' | |
3680 | c WRITE(ifps,'(a)') '%%BeginSetup' | |
3681 | c WRITE(ifps,'(a)') '%%IncludeFeature: *PageSize A4' | |
3682 | c WRITE(ifps,'(a)') '%%EndSetup' | |
3683 | c WRITE(ifps,'(a)') '%%EndComments' | |
3684 | c WRITE(ifps,*) '/l {lineto} bind def' | |
3685 | c WRITE(ifps,*) '/rl {rlineto} bind def' | |
3686 | c WRITE(ifps,*) '/m {moveto} bind def' | |
3687 | c WRITE(ifps,*) '/rm {rmoveto} bind def' | |
3688 | c WRITE(ifps,*) '/s {stroke} bind def' | |
3689 | c WRITE(ifps,*) '/gr {grestore} bind def' | |
3690 | c WRITE(ifps,*) '/gs {gsave} bind def' | |
3691 | c WRITE(ifps,*) '/cp {closepath} bind def' | |
3692 | c WRITE(ifps,*) '/tr {translate} bind def' | |
3693 | c WRITE(ifps,*) '/sc {scale} bind def' | |
3694 | c WRITE(ifps,*) '/sd {setdash} bind def' | |
3695 | c WRITE(ifps,*) '/sdo {[.01 .05] 0 sd} bind def' | |
3696 | c WRITE(ifps,*) '/sdf {[1 .0] 0 sd} bind def' | |
3697 | c WRITE(ifps,*) '/n {newpath} bind def' | |
3698 | c WRITE(ifps,*) '/slw {setlinewidth } bind def' | |
3699 | c write(ifps,*) '/srgb {setrgbcolor} bind def' | |
3700 | c write(ifps,*) '/lgrey { 0 0.95 0.95 srgb} bind def' | |
3701 | c write(ifps,*) '/black { 0 0 0 srgb} bind def' | |
3702 | c write(ifps,*) '/red { 1 0 0 srgb} bind def ' | |
3703 | c write(ifps,*) '/green { 0 1 0 srgb} bind def ' | |
3704 | c write(ifps,*) '/blue { 0 0 1 srgb} bind def ' | |
3705 | c write(ifps,*) '/yellow { 1 0.5 0 srgb} bind def ' | |
3706 | c write(ifps,*) '/turquoise { 0 1 1 srgb} bind def ' | |
3707 | c write(ifps,*) '/purple { 1 0 1 srgb} bind def ' | |
3708 | cc write(ifps,*) '/ { srgb} bind def ' | |
3709 | cc write(ifps,*) '/ { srgb} bind def ' | |
3710 | c write(ifps,*) '/ef {eofill} bind def' | |
3711 | c WRITE(ifps,'(a)') '%%EndProlog' | |
3712 | c WRITE(ifps,*) 'gsave' | |
3713 | c WRITE(ifps,*) '/Helvetica findfont 10 scalefont setfont' | |
3714 | c color(9)='lgrey ' | |
3715 | c color(1)='black ' | |
3716 | c color(2)='red ' | |
3717 | c color(3)='green ' | |
3718 | c color(4)='blue ' | |
3719 | c color(7)='yellow ' | |
3720 | c color(5)='turquoise ' | |
3721 | c color(6)='purple ' | |
3722 | c np=0 | |
3723 | c elseif(iii.eq.1)then | |
3724 | c np=np+1 | |
3725 | c write(ifps,'(a,i4)') '%%Page: number ',np | |
3726 | c write(ifps,'(a)') 'gsave' | |
3727 | c WRITE(ifps,*) '100 700 tr' | |
3728 | c scale=0.125 | |
3729 | c WRITE(ifps,*) 1./scale,1./scale,' sc' | |
3730 | c WRITE(ifps,*) scale/2.,' slw' | |
3731 | c WRITE(ifps,*) '/Helvetica findfont ',15.*scale | |
3732 | c & ,' scalefont setfont' | |
3733 | c write(ifps,*) color(1),' n ',smin,xmin,' m ( tau:',tau,') show ' | |
3734 | c | |
3735 | c WRITE(ifps,*) '/Helvetica findfont ',5.*scale | |
3736 | c & ,' scalefont setfont' | |
3737 | c | |
3738 | c | |
3739 | c yb=-2. | |
3740 | c dy=4./12. | |
3741 | c yb=yb-dy/2 | |
3742 | c do iyb=0,11 | |
3743 | c yb=yb+dy | |
3744 | c WRITE(ifps,*) 'gsave' | |
3745 | c WRITE(ifps,*) (xmax-xmin)*1.1*float(int(iyb/4)) | |
3746 | c & ,-(xmax-xmin)*1.1*mod(iyb,4),' tr' | |
3747 | c write(ifps,*) ' n ',smin,xmin,' m ',smax,xmin,' l ' | |
3748 | c & ,smax,xmax,' l ',smin,xmax,' l cp s ' | |
3749 | cc.......particles in layer iyb............. | |
3750 | c do i=1,nptl | |
3751 | c if(ii.gt.0)then | |
3752 | c write(ifps,*) color(mod(i,5)+2) | |
3753 | c & ,' n ',u,x-r,' m ',o,x-r,' l ' | |
3754 | c & ,o,x+r,' l ',u,x+r,' l cp s ' | |
3755 | c write(ifps,*) ' n ',u,x-r,' m (',i,ior,') show ' | |
3756 | c else | |
3757 | c write(ifps,*) ' n ',s,x,r,0,360,' arc ',color(ic),' s ' | |
3758 | c write(ifps,*) ' n ',s-r,x,' m (',i,io,') show ' | |
3759 | c endif | |
3760 | c 10 enddo | |
3761 | c write(ifps,*) color(1),' n ',smin,xmin,' m (',yb,') show ' | |
3762 | c WRITE(ifps,*) 'grestore' | |
3763 | c enddo !yb bin | |
3764 | c write(ifps,'(a)') 'grestore' | |
3765 | c write(ifps,*) 'showpage' | |
3766 | c elseif(iii.eq.2)then | |
3767 | c write(ifps,*) 'gr' | |
3768 | c | |
3769 | c write(ifps,'(a)') '%%Trailer' | |
3770 | c write(ifps,'(a,i4)') '%%Pages: ',np | |
3771 | c write(ifps,'(a)') '%%EOF' | |
3772 | c close(unit=ifps) | |
3773 | c endif | |
3774 | c | |
3775 | c return | |
3776 | c end | |
3777 | c | |
3778 | ||
3779 | c---------------------------------------------------------------------------------------- | |
3780 | subroutine xtauev(iii) | |
3781 | c---------------------------------------------------------------------------------------- | |
3782 | end | |
3783 | c---------------------------------------------------------------------------------------- | |
3784 | subroutine wimi | |
3785 | c---------------------------------------------------------------------------------------- | |
3786 | end | |
3787 | c---------------------------------------------------------------------------------------- | |
3788 | subroutine wimino | |
3789 | c---------------------------------------------------------------------------------------- | |
3790 | end | |
3791 | c---------------------------------------------------------------------------------------- | |
3792 | subroutine xspace(iii) | |
3793 | c---------------------------------------------------------------------------------------- | |
3794 | end | |
3795 | c---------------------------------------------------------------------------------------- | |
3796 | subroutine wclu | |
3797 | c---------------------------------------------------------------------------------------- | |
3798 | end | |
3799 | c---------------------------------------------------------------------------------------- | |
3800 | subroutine wclufi | |
3801 | c---------------------------------------------------------------------------------------- | |
3802 | end | |
3803 | c---------------------------------------------------------------------------------------- | |
3804 | subroutine wtime(iii) | |
3805 | c---------------------------------------------------------------------------------------- | |
3806 | end |