]> git.uio.no Git - u/mrichter/AliRoot.git/blob - EPOS/epos167/epos-int-168.f
Fixes needed by gfortran, it doesn't perform lazy evaluation (Piotr)
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-int-168.f
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