Removing leftover return

[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-int-168.f

c-----------------------------------------------------------------------
      subroutine bjinta(ier)
c-----------------------------------------------------------------------
c  fin. state interactions and decays
c-----------------------------------------------------------------------
      include 'epos.inc'
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctimel/ntc
      common/col3/ncol,kolpt
      double precision ttaun,ttau0,rcproj,rctarg
      common/cttaun/ttaun /cttau0/ttau0 /geom1/rcproj,rctarg
      logical go

      call utpri('bjinta',ish,ishini,4)
     
      ier=0
      
      if(ncol.eq.0.and.iappl.eq.2)goto1000
      if(nevt.ne.1.or.ifrade.eq.0)goto1000

      if(iappl.eq.4)then
ctp060829        iacn=0
        goto5000
      endif

      if(iappl.eq.1)then
        tauxx=0.7+0.94*max(radnuc(maproj),radnuc(matarg))/(0.5*engy)
      else
        tauxx=0.
      endif
      tauzz=max(taumin,tauxx)
      !print*,'====',taumin,tauxx,tauzz
      ttaus=dble(tauzz)
      ttau0=dsqrt(rcproj*rctarg)
      call jtauin   ! initialize hyperbola

      if(iappl.ne.1)goto 5000


c     no-secondary-interactions or parton-ladder-fusion
c     -------------------------------------------------
      if(iorsce.eq.0.and.iorsdf.eq.0.and.iorshh.eq.0
     &      .or.iorsdf.eq.3)then
ctp060829        iacn=0
        if(iorsdf.eq.3)then
c if nptl already very big, clean up useless particles in cptl list.
c (do not use it when gakstr() is called (some information lost)
          if(nclean.gt.0.and.nptl.gt.mxptl/5)then
            nptli=maproj+matarg+1
            do iii=nptli,nptl
              go=.true.
              if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
              if(go.and.mod(istptl(iii),10).ne.0)istptl(iii)=99
            enddo
            nptl0=nptl
            call utclea(nptli,nptl0)
          endif
          nptlbpo=nptl
          call jintpo  !parton-ladder-fusion
         if(ish.ge.2)call alist('parton-ladder-fusion&',nptlbpo+1,nptl)
        endif
        goto5000
      else
        stop'bjinta: not supported any more (310305).     '
      endif

5000  continue

      nptlbd=nptl
      
      call xSpaceTime

      if(ifrade.eq.0)goto779  !skip decay
      if(idecay.eq.0)goto779  !skip decay


      if(ish.ge.2)call alist('final decay&',0,0)
      if(iappl.eq.7)then
        nptli=1
      else
        nptli=maproj+matarg+1
      endif
      np1=nptli
41    np2=nptl
      nptli=np1
      ip=np1-1
      do while (ip.lt.np2)
      ip=ip+1
      if(istptl(ip).eq.0)then
      call hdecas(ip,iret)
      if(iret.eq.1)goto 1001
      if(iret.eq.-1)goto 42
c remove useless particles if not enough space
      if(nclean.gt.0.and.nptl.gt.mxptl/2)then
        nnnpt=0
        do iii=nptli,ip
          go=.true.
          if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
          if(go.and.mod(istptl(iii),10).ne.0)then
            istptl(iii)=99
            nnnpt=nnnpt+1
          endif
        enddo
        if(nnnpt.gt.mxptl-nptl)then
          nptl0=nptl
          call utclea(nptli,nptl0)
          np2=np2-nnnpt
          ip=ip-nnnpt
          nptli=ip
        endif
      endif
      endif
42    continue
      enddo
      nptli=max(nptli,np1)
      np1=np2+1
      if(np1.le.nptl)then
      if(ish.ge.2)then
      if(ish.ge.3)call alist('partial list',0,0)
      do 6 ip=np1,nptl
        call alist('&',ip,ip)
6     continue
      endif
      goto 41
      endif
  779 continue

      if(ish.ge.2)call alist('complete list&',1,nptl)
      
c     on shell check
c     --------------
c      if(iappl.eq.1)call jresc
      
1000  continue
      call utprix('bjinta',ish,ishini,4)
      return

1001  continue
      ier=1
      goto1000

      end

cc----------------------------------------------------------------------
c      subroutine jintcs(i,j,ecm,bij,nq,jc,ics)
cc----------------------------------------------------------------------
cc compare hadron distance with energy dependent cross section
cc data taken from particle data group, durham and juelich
cc input:
cc   i,j: particle indices
cc   ecm: center-of-mass energy  
cc   bij: impact parameter
cc   nq: net quark number of fused object
cc   jc: jc of fused object
cc output: 
cc   ics=0 if distance larger than sqrt(sig(E_CMS)/pi)
cc   ics=1 else
cc The data are from HEPDATA, 
cc the formulas from Rev. Particle Properties 1995
cc----------------------------------------------------------------------
c      include 'epos.inc'
c      integer jci(nflav,2),jcj(nflav,2),jc(nflav,2),kc(nflav)
c     *,kci(nflav),kcj(nflav)
c      common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
c     *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
c      common/cdfptl/idfptl(mxptl)
c      parameter(npp=249,napp=205,npn=411,napn=31,npip=441)
c      parameter(npim=578,nkmp=299,nkmn=41,nkpp=172,nkpn=91)
cc      parameter(npim=578,nkmp=299,nkmn=41,nkpp=172,nkpn=91,nlp=35)
c      parameter(npi1=12,npi2=12,npi3=18,npi4=21,npi5=9)
cc      real ppecm(npp)
c      real ppbmx(npp)
c      real appecm(napp),appbmx(napp)
c      real pnecm(npn),pnbmx(npn)
c      real apnecm(napn),apnbmx(napn)
c      real pipecm(npip),pipbmx(npip)
c      real pimecm(npim),pimbmx(npim)
c      real kmpecm(nkmp),kmpbmx(nkmp)
c      real kmnecm(nkmn),kmnbmx(nkmn)
c      real kppecm(nkpp),kppbmx(nkpp)
c      real kpnecm(nkpn),kpnbmx(nkpn)
cc      real lpecm(nlp),lpbmx(nlp)
c      real pi1ecm(npi1),pi1bmx(npi1)
c      real pi2ecm(npi2),pi2bmx(npi2)
c      real pi3ecm(npi3),pi3bmx(npi3)
c      real pi4ecm(npi4),pi4bmx(npi4)
c      real pi5ecm(npi5),pi5bmx(npi5)
c
cc      data ppecm/
cc     *    1.8812,   1.8855,   1.8910,   1.8963,   1.9073
cc     *,   1.9108,   1.9145,   1.9224,   1.9244,   1.9352
cc     *,   1.9466,   1.9468,   1.9542,   1.9592,   1.9636
cc     *,   1.9772,   1.9860,   1.9945,   2.0032,   2.0052
cc     *,   2.0070,   2.0272,   2.0275,   2.0302,   2.0333
cc     *,   2.0402,   2.0427,   2.0586,   2.0608,   2.0692
cc     *,   2.0702,   2.0708,   2.0715,   2.0718,   2.0751
cc     *,   2.0797,   2.0813,   2.0829,   2.0843,   2.0846
cc     *,   2.0935,   2.1062,   2.1113,   2.1123,   2.1170
cc     *,   2.1173,   2.1184,   2.1268,   2.1289,   2.1357
cc     *,   2.1467,   2.1511,   2.1522,   2.1553,   2.1618
cc     *,   2.1639,   2.1726,   2.1771,   2.1795,   2.1799
cc     *,   2.1802,   2.1813,   2.1868,   2.2152,   2.2184
cc     *,   2.2212,   2.2254,   2.2395,   2.2405,   2.2532
cc     *,   2.2606,   2.2613,   2.2861,   2.2889,   2.2914
cc     *,   2.2988,   2.3101,   2.3136,   2.3228,   2.3348
cc     *,   2.3490,   2.3525,   2.3529,   2.3800,   2.3842
cc     *,   2.3912,   2.4088,   2.4130,   2.4193,   2.4298
cc     *,   2.4315,   2.4392,   2.4472,   2.4573,   2.5034
cc     *,   2.5131,   2.5268,   2.5743,   2.5848,   2.5916
cc     *,   2.6327,   2.6620,   2.6700,   2.6984,   2.7080
cc     *,   2.7205,   2.7534,   2.7573,   2.7651,   2.7670
cc     *,   2.7844,   2.8024,   2.8092,   2.8127,   2.8533
cc     *,   2.8556,   2.8638,   2.9079,   2.9395,   2.9500
cc     *,   2.9776,   2.9961,   3.0495,   3.0769,   3.0879
cc     *,   3.1358,   3.1547,   3.2251,   3.2371,   3.3020
cc     *,   3.3527,   3.3620,   3.4221,   3.4967,   3.5019
cc     *,   3.5035,   3.5814,   3.5829,   3.6266,   3.8549
cc     *,   3.8742,   4.0503,   4.0663,   4.0698,   4.0732
cc     *,   4.0778,   4.1074,   4.1300,   4.4976,   4.5183
cc     *,   4.5389,   4.5410,   4.5615,   4.6808,   4.9146
cc     *,   4.9336,   5.0088,   5.2993,   5.3011,   5.4731
cc     *,   5.6083,   5.6416,   5.6465,   5.9171,   5.9502
cc     *,   5.9644,   6.1697,   6.1803,   6.2706,   6.3034
cc     *,   6.3390,   6.5627,   6.7040,   6.8424,   6.9105
cc     *,   6.9780,   7.1111,   7.1662,   7.6202,   7.8624
cc     *,   8.2124,   8.7647,   9.0282,   9.2843,   9.5825
cc     *,   9.7763,   9.9851,  10.2447,  10.6927,  10.8926
cc     *,  11.4549,  11.5365,  11.7779,  11.8519,  13.6241
cc     *,  13.6883,  13.7611,  13.8968,  15.0628,  15.1868
cc     *,  16.6595,  16.8275,  17.9077,  18.0121,  18.1677
cc     *,  19.2213,  19.4156,  19.6556,  19.7002,  21.2604
cc     *,  22.9574,  23.3624,  23.4057,  23.4965,  23.4967
cc     *,  23.5964,  23.7605,  23.8787,  24.1521,  25.2904
cc     *,  26.3796,  27.5960,  30.5240,  30.5954,  30.5957
cc     *,  30.6555,  30.6954,  30.7954,  35.1947,  44.6933
cc     *,  44.6937,  44.6938,  44.7706,  44.8228,  44.8933
cc     *,  45.1933,  52.6798,  52.7090,  52.7921,  52.7921
cc     *,  52.7927,  52.8927,  53.1921,  62.2907,  62.3914
cc     *,  62.4907,  62.6907,  62.6913,  62.7906
cc     */
c      data ppbmx/
c     *    3.1615,   2.2212,   1.7113,   1.4927,   1.1631
c     *,   1.0911,   1.0388,    .9525,    .9390,    .8885
c     *,    .8019,    .8956,    .9115,    .8593,    .8840
c     *,    .9062,    .8444,    .8444,    .8482,    .8713
c     *,    .8575,    .8795,    .8795,    .8759,    .8630
c     *,    .8822,    .8593,    .8813,    .9036,    .9150
c     *,    .8740,    .9219,    .8740,    .9253,    .8722
c     *,    .8795,    .9288,    .8704,    .9398,    .9271
c     *,    .9407,    .9373,    .9756,    .9926,   1.0357
c     *,   1.0037,   1.0408,    .9739,   1.0108,   1.0479
c     *,   1.0645,   1.0705,   1.1113,   1.0794,   1.0955
c     *,   1.1085,   1.1256,   1.1535,   1.1731,   1.1424
c     *,   1.0794,   1.1480,   1.1590,   1.1888,   1.2218
c     *,   1.2164,   1.2087,   1.2296,   1.2231,   1.2335
c     *,   1.2322,   1.2309,   1.2114,   1.2296,   1.2293
c     *,   1.2489,   1.2303,   1.2322,   1.2322,   1.2127
c     *,   1.2192,   1.2295,   1.2399,   1.2290,   1.2374
c     *,   1.2205,   1.2278,   1.2284,   1.2061,   1.2296
c     *,   1.2296,   1.2540,   1.2008,   1.2260,   1.2229
c     *,   1.2257,   1.2188,   1.2118,   1.2078,   1.1982
c     *,   1.2039,   1.2012,   1.1991,   1.1808,   1.1969
c     *,   1.1959,   1.1922,   1.1902,   1.1897,   1.1878
c     *,   1.1888,   1.1860,   1.1856,   1.1850,   1.2244
c     *,   1.1782,   1.1790,   1.1718,   1.1696,   1.1726
c     *,   1.1576,   1.1656,   1.1606,   1.1603,   1.1581
c     *,   1.1590,   1.1530,   1.1576,   1.1487,   1.1476
c     *,   1.1447,   1.1794,   1.1448,   1.1507,   1.1507
c     *,   1.1407,   1.1403,   1.1507,   1.1581,   1.1645
c     *,   1.1616,   1.1507,   1.1332,   1.1294,   1.1284
c     *,   1.1227,   1.1284,   1.1298,   1.1261,   1.1199
c     *,   1.1365,   1.1438,   1.1284,   1.1424,   1.1230
c     *,   1.1218,   1.1142,   1.1156,   1.1202,   1.1198
c     *,   1.1099,   1.1099,   1.1175,   1.1241,   1.1168
c     *,   1.1099,   1.1128,   1.1241,   1.1103,   1.1149
c     *,   1.1155,   1.1083,   1.1197,   1.1127,   1.1185
c     *,   1.1113,   1.1128,   1.1113,   1.1083,   1.1056
c     *,   1.1067,   1.1070,   1.1059,   1.1063,   1.1070
c     *,   1.1028,   1.0979,   1.1060,   1.1062,   1.1774
c     *,   1.0805,   1.1039,   1.0940,   1.0955,   1.0940
c     *,   1.0955,   1.1082,   1.1128,   1.1082,   1.0852
c     *,   1.1003,   1.1097,   1.1118,   1.0911,   1.1227
c     *,   1.1070,   1.1206,   1.1142,   1.1213,   1.1174
c     *,   1.1199,   1.1142,   1.1185,   1.1180,   1.1113
c     *,   1.1099,   1.1297,   1.1142,   1.1226,   1.1240
c     *,   1.1251,   1.1368,   1.1403,   1.1319,   1.1296
c     *,   1.1333,   1.1303,   1.1284,   1.1343,   1.1516
c     *,   1.1521,   1.1549,   1.1641,   1.1631,   1.1562
c     *,   1.1631,   1.1697,   1.1726,   1.1756,   1.1659
c     *,   1.1660,   1.1617,   1.1686,   1.1774,   1.1713
c     *,   1.1744,   1.1810,   1.1694,   1.1848
c     */
c      data appecm/
c     *    1.9002,   1.9050,   1.9072,   1.9078,   1.9091
c     *,   1.9129,   1.9157,   1.9162,   1.9174,   1.9176
c     *,   1.9180,   1.9195,   1.9201,   1.9224,   1.9226
c     *,   1.9246,   1.9252,   1.9255,   1.9257,   1.9271
c     *,   1.9282,   1.9293,   1.9301,   1.9310,   1.9319
c     *,   1.9328,   1.9334,   1.9345,   1.9359,   1.9370
c     *,   1.9372,   1.9384,   1.9393,   1.9398,   1.9407
c     *,   1.9426,   1.9430,   1.9433,   1.9452,   1.9454
c     *,   1.9473,   1.9485,   1.9495,   1.9500,   1.9510
c     *,   1.9515,   1.9547,   1.9559,   1.9562,   1.9579
c     *,   1.9610,   1.9615,   1.9644,   1.9680,   1.9718
c     *,   1.9755,   1.9788,   1.9829,   1.9871,   1.9911
c     *,   1.9954,   1.9994,   2.0813,   2.0979,   2.1146
c     *,   2.1180,   2.1316,   2.1487,   2.1660,   2.1834
c     *,   2.1868,   2.1938,   2.1991,   2.2008,   2.2184
c     *,   2.2226,   2.2359,   2.2500,   2.2606,   2.2712
c     *,   2.2889,   2.2995,   2.3066,   2.3243,   2.3419
c     *,   2.3490,   2.3511,   2.3596,   2.3701,   2.3772
c     *,   2.3860,   2.3877,   2.3947,   2.4035,   2.4123
c     *,   2.4298,   2.4472,   2.4629,   2.4820,   2.4993
c     *,   2.5165,   2.5268,   2.5337,   2.5508,   2.5678
c     *,   2.5848,   2.5950,   2.6017,   2.6186,   2.6353
c     *,   2.6377,   2.6520,   2.6654,   2.6687,   2.6852
c     *,   2.7017,   2.7182,   2.7345,   2.7508,   2.7670
c     *,   2.7832,   2.7896,   2.7992,   2.8152,   2.8312
c     *,   2.8439,   2.8470,   2.8565,   2.8628,   2.9377
c     *,   2.9561,   2.9745,   3.0351,   3.0769,   3.0828
c     *,   3.1648,   3.2507,   3.2788,   3.3620,   3.4568
c     *,   3.5492,   3.6266,   3.7893,   3.8501,   3.8597
c     *,   3.8742,   3.9455,   4.1074,   4.3499,   4.3926
c     *,   4.5389,   4.5799,   4.6808,   4.7991,   4.9336
c     *,   4.9901,   5.1742,   5.2993,   5.3520,   5.4731
c     *,   5.6416,   5.6911,   5.8534,   5.9644,   6.0113
c     *,   6.2706,   6.3153,   6.7040,   6.9780,   7.3061
c     *,   7.6202,   7.7422,   7.8624,   7.8743,   8.0393
c     *,   8.2124,   8.4930,   8.7647,   9.0282,   9.2843
c     *,   9.5335,   9.7763,  11.5365,  13.7611,  13.8630
c     *,  15.0628,  16.8275,  17.9077,  19.4156,  21.2604
c     *,  22.9574,  30.4098,  30.5943,  30.6861,  52.5843
c     *,  52.7979,  52.7979,  62.2853,  62.4957,  62.6905
c     *, 539.9198, 546.9191, 899.8658, 900.0000,1803.0007
c     */
c      data appbmx/
c     *    2.7553,   2.6857,   2.6732,   2.6517,   2.6306
c     *,   2.5376,   2.4398,   2.4779,   2.5212,   2.5143
c     *,   2.4573,   2.4799,   2.4495,   2.4286,   2.4521
c     *,   2.3683,   2.4463,   2.4240,   2.4489,   2.3743
c     *,   2.4162,   2.3561,   2.3903,   2.3796,   2.3221
c     *,   2.3689,   2.4181,   2.3221,   2.3823,   2.2987
c     *,   2.3337,   2.3568,   2.3487,   2.2911,   2.3344
c     *,   2.2701,   2.3097,   2.3207,   2.2631,   2.2603
c     *,   2.2666,   2.2412,   2.2319,   2.3097,   2.2327
c     *,   2.2104,   2.2104,   2.2645,   2.2369,   2.1880
c     *,   2.1756,   2.2162,   2.1491,   2.1484,   2.1365
c     *,   2.1320,   2.1027,   2.1065,   2.0913,   2.0776
c     *,   2.0668,   2.0883,   1.9333,   1.9098,   1.8851
c     *,   1.8932,   1.8721,   1.8623,   1.8520,   1.8409
c     *,   1.8678,   1.8575,   1.8325,   1.8429,   1.8088
c     *,   1.7698,   1.7941,   1.7864,   1.7814,   1.7736
c     *,   1.7645,   1.7623,   1.7576,   1.7523,   1.7445
c     *,   1.7419,   1.7019,   1.7342,   1.7311,   1.7271
c     *,   1.7210,   1.7182,   1.7161,   1.7119,   1.7054
c     *,   1.6947,   1.6816,   1.6780,   1.6678,   1.6599
c     *,   1.6509,   1.6328,   1.6449,   1.6396,   1.6319
c     *,   1.6270,   1.5761,   1.6043,   1.6120,   1.6069
c     *,   1.5978,   1.6018,   1.5891,   1.5948,   1.5905
c     *,   1.5849,   1.5778,   1.5736,   1.5679,   1.5785
c     *,   1.5583,   1.5476,   1.5519,   1.5467,   1.5416
c     *,   1.5233,   1.5368,   1.5492,   1.5313,   1.4895
c     *,   1.5579,   1.5107,   1.4680,   1.4725,   1.4586
c     *,   1.3889,   1.4799,   1.4472,   1.4371,   1.3576
c     *,   1.4228,   1.3922,   1.3762,   1.3854,   1.3669
c     *,   1.4161,   1.3623,   1.3638,   1.3530,   1.3762
c     *,   1.3273,   1.2866,   1.2989,   1.3159,   1.2828
c     *,   1.3086,   1.2797,   1.2704,   1.2940,   1.2741
c     *,   1.2514,   1.2540,   1.2803,   1.2653,   1.2603
c     *,   1.2386,   1.2192,   1.2283,   1.2231,   1.2140
c     *,   1.2048,   1.2114,   1.2099,   1.2052,   1.2048
c     *,   1.1984,   1.1928,   1.1901,   1.1902,   1.1888
c     *,   1.1848,   1.1769,   1.1706,   1.1579,   1.1608
c     *,   1.1521,   1.1534,   1.1520,   1.1495,   1.1549
c     *,   1.1550,   1.1580,   1.1672,   1.1562,   1.1743
c     *,   1.1859,   1.1950,   1.1851,   1.1821,   1.1987
c     *,   1.4740,   1.4037,   1.4483,   1.4417,   1.5149
c     */
c      data pnecm/
c     *    1.87867,   1.87869,   1.87872,   1.87875,   1.87878
c     *,   1.87881,   1.87884,   1.87887,   1.87890,   1.87893
c     *,   1.87896,   1.87899,   1.87902,   1.87906,   1.87909
c     *,   1.87913,   1.87916,   1.87920,   1.87923,   1.87927
c     *,   1.87931,   1.87934,   1.87938,   1.87942,   1.87946
c     *,   1.87950,   1.87954,   1.87958,   1.87962,   1.87966
c     *,   1.87970,   1.87975,   1.87979,   1.87983,   1.87988
c     *,   1.87992,   1.87997,   1.88001,   1.88006,   1.88011
c     *,   1.88015,   1.88020,   1.88025,   1.88030,   1.88035
c     *,   1.88040,   1.88045,   1.88050,   1.88055,   1.88061
c     *,   1.88066,   1.88071,   1.88077,   1.88082,   1.88087
c     *,   1.88093,   1.88099,   1.88104,   1.88110,   1.88116
c     *,   1.88121,   1.88127,   1.88133,   1.88139,   1.88145
c     *,   1.88151,   1.88157,   1.88163,   1.88170,   1.88176
c     *,   1.88182,   1.88189,   1.88195,   1.88202,   1.88208
c     *,   1.88215,   1.88221,   1.88228,   1.88235,   1.88241
c     *,   1.88248,   1.88255,   1.88262,   1.88269,   1.88276
c     *,   1.88283,   1.88290,   1.88297,   1.88305,   1.88312
c     *,   1.88319,   1.88327,   1.88334,   1.88342,   1.88349
c     *,   1.88357,   1.88364,   1.88372,   1.88380,   1.88388
c     *,   1.88396,   1.88403,   1.88411,   1.88419,   1.88428
c     *,   1.88436,   1.88444,   1.88452,   1.88460,   1.88469
c     *,   1.88477,   1.88485,   1.88494,   1.88502,   1.88511
c     *,   1.88519,   1.88528,   1.88537,   1.88546,   1.88554
c     *,   1.88563,   1.88572,   1.88581,   1.88590,   1.88599
c     *,   1.88608,   1.88618,   1.88627,   1.88636,   1.88645
c     *,   1.88655,   1.88664,   1.88674,   1.88683,   1.88693
c     *,   1.88702,   1.88712,   1.88722,   1.88731,   1.88741
c     *,   1.88751,   1.88761,   1.88771,   1.88781,   1.88791
c     *,   1.88801,   1.88811,   1.88822,   1.88832,   1.88842
c     *,   1.88852,   1.88863,   1.88873,   1.88884,   1.88894
c     *,   1.88905,   1.88916,   1.88926,   1.88937,   1.88948
c     *,   1.88959,   1.88970,   1.88980,   1.88991,   1.89003
c     *,   1.89014,   1.89025,   1.89036,   1.89047,   1.89058
c     *,   1.89070,   1.89081,   1.89093,   1.89104,   1.89116
c     *,   1.89127,   1.89139,   1.89150,   1.89162,   1.89174
c     *,   1.89186,   1.89198,   1.89209,   1.89221,   1.89233
c     *,   1.89245,   1.89258,   1.89270,   1.89282,   1.89294
c     *,   1.89306,   1.89319,   1.89331,   1.89344,   1.89356
c     *,   1.89369,   1.89381,   1.89394,   1.89406,   1.89419
c     *,   1.89432,   1.89445,   1.89458,   1.89496,   1.89549
c     *,   1.89602,   1.89615,   1.89656,   1.89697,   1.89724
c     *,   1.89766,   1.89780,   1.89850,   1.89893,   1.89922
c     *,   1.89995,   1.90068,   1.90098,   1.90143,   1.90174
c     *,   1.90235,   1.90250,   1.90312,   1.90407,   1.90503
c     *,   1.90519,   1.90616,   1.90732,   1.90749,   1.90833
c     *,   1.90936,   1.91216,   1.91379,   1.91545,   1.91714
c     *,   1.91905,   1.91924,   1.92100,   1.92160,   1.92179
c     *,   1.92259,   1.92319,   1.92421,   1.92502,   1.92543
c     *,   1.92605,   1.92646,   1.92792,   1.93046,   1.93241
c     *,   1.93306,   1.93570,   1.93704,   1.93953,   1.94021
c     *,   1.94183,   1.94699,   1.94723,   1.95206,   1.95329
c     *,   1.95452,   1.95651,   1.96029,   1.96080,   1.97927
c     *,   1.98533,   1.99806,   2.00884,   2.01269,   2.02779
c     *,   2.02963,   2.04329,   2.04643,   2.05915,   2.05947
c     *,   2.05979,   2.07207,   2.07337,   2.07533,   2.09178
c     *,   2.10847,   2.11385,   2.12061,   2.12536,   2.14243
c     *,   2.15308,   2.15343,   2.15964,   2.17072,   2.17698
c     *,   2.18185,   2.19442,   2.21193,   2.21369,   2.22353
c     *,   2.22846,   2.22952,   2.24715,   2.25597,   2.28851
c     *,   2.29134,   2.29382,   2.31293,   2.31328,   2.32673
c     *,   2.34935,   2.35501,   2.36207,   2.38251,   2.39729
c     *,   2.41135,   2.41591,   2.42817,   2.44006,   2.45995
c     *,   2.48220,   2.49952,   2.50643,   2.52951,   2.56750
c     *,   2.58756,   2.63546,   2.64719,   2.66487,   2.67285
c     *,   2.71088,   2.72336,   2.75178,   2.75634,   2.76802
c     *,   2.76867,   2.76996,   2.78741,   2.80542,   2.81567
c     *,   2.85639,   2.85860,   2.86681,   2.91102,   2.94265
c     *,   2.96543,   3.05273,   3.08729,   3.09115,   3.15805
c     *,   3.16821,   3.22857,   3.24053,   3.29551,   3.35628
c     *,   3.42580,   3.50727,   3.55297,   3.58520,   3.58675
c     *,   3.59580,   3.63049,   3.72425,   3.75633,   4.05460
c     *,   4.06219,   4.07412,   4.11175,   4.37426,   4.52314
c     *,   4.54378,   4.74537,   4.93885,   5.12513,   5.30495
c     *,   5.40999,   5.47892,   5.61424,   5.64757,   5.93444
c     *,   5.97071,   6.27732,   6.42517,   6.51227,   6.56970
c     *,   6.71114,   6.87703,   6.98545,   7.18446,   7.23904
c     *,   7.24942,   7.62830,   8.10603,   8.22114,   8.54948
c     *,   8.77406,   9.29418,   9.78673,  10.15712,  10.25566
c     *,  10.70407,  11.13446,  11.54882,  12.33587,  13.77577
c     *,  15.07880,  15.74960,  16.84544,  17.92675,  18.18703
c     *,  18.44365,  19.43624,  20.14863,  21.28302,  22.69375
c     *,  22.98187
c     */
c      data pnbmx/
c     *   10.7613,  10.6195,  10.5290,  10.4010,  10.2901
c     *,  10.1970,  10.1029,  10.0088,   9.9009,   9.7973
c     *,   9.7125,   9.6198,   9.5613,   9.4281,   9.3628
c     *,   9.2801,   9.2068,   9.1103,   9.0239,   8.9609
c     *,   8.8562,   8.8011,   8.7574,   8.6332,   8.5487
c     *,   8.4788,   8.4157,   8.3606,   8.2746,   8.2604
c     *,   8.1583,   8.0778,   8.0127,   7.9567,   7.8691
c     *,   7.8176,   7.8190,   7.7028,   7.6634,   7.5956
c     *,   7.5808,   7.4956,   7.4144,   7.3600,   7.3071
c     *,   7.2947,   7.1947,   7.1389,   7.1333,   7.0404
c     *,   6.9856,   6.9468,   6.9171,   6.8369,   6.7869
c     *,   6.7829,   6.6804,   6.6483,   6.5982,   6.5464
c     *,   6.5060,   6.4541,   6.4539,   6.3684,   6.3219
c     *,   6.2777,   6.2832,   6.1881,   6.1494,   6.1048
c     *,   6.0642,   6.0207,   6.0255,   5.9756,   5.9161
c     *,   5.8717,   5.8276,   5.7878,   5.7481,   5.7319
c     *,   5.6977,   5.6492,   5.6142,   5.5796,   5.5428
c     *,   5.4939,   5.4659,   5.4549,   5.3933,   5.3601
c     *,   5.3311,   5.3037,   5.2760,   5.2429,   5.2381
c     *,   5.2384,   5.1470,   5.1112,   5.0812,   5.0483
c     *,   5.0350,   5.0330,   4.9640,   4.9410,   4.9050
c     *,   4.8714,   4.8408,   4.8482,   4.8518,   4.7611
c     *,   4.7157,   4.7049,   4.6834,   4.6939,   4.6206
c     *,   4.5931,   4.5723,   4.5498,   4.5210,   4.5337
c     *,   4.4726,   4.4898,   4.4251,   4.4013,   4.3784
c     *,   4.3514,   4.3264,   4.3530,   4.2846,   4.2567
c     *,   4.2339,   4.2171,   4.1919,   4.1757,   4.1503
c     *,   4.1329,   4.0924,   4.0797,   4.0624,   4.0446
c     *,   4.0210,   4.0010,   3.9812,   3.9627,   3.9402
c     *,   3.9221,   3.9029,   3.8925,   3.8680,   3.8439
c     *,   3.8244,   3.8094,   3.7816,   3.7719,   3.7396
c     *,   3.7271,   3.7011,   3.6849,   3.6804,   3.6640
c     *,   3.6496,   3.5954,   3.6042,   3.5697,   3.5552
c     *,   3.5463,   3.5265,   3.5212,   3.4992,   3.1432
c     *,   3.4630,   3.4483,   3.4484,   3.4240,   3.4127
c     *,   3.3595,   3.3767,   3.3582,   3.3059,   3.3258
c     *,   3.3182,   3.2582,   3.2740,   3.2621,   3.1984
c     *,   3.2408,   3.2043,   3.2087,   3.1652,   3.1720
c     *,   3.1476,   3.0754,   3.1236,   3.1494,   3.0653
c     *,   3.0757,   3.1030,   3.0692,   3.0466,   3.0392
c     *,   3.0346,   2.9933,   3.0146,   2.8774,   2.7972
c     *,   2.6463,   2.8068,   2.6869,   2.6643,   2.6673
c     *,   2.6684,   2.6409,   2.5520,   2.4341,   2.4978
c     *,   2.4534,   2.4417,   2.4293,   2.3283,   2.4978
c     *,   2.3042,   2.3111,   2.1996,   2.1996,   2.1960
c     *,   2.1469,   2.1223,   2.0791,   2.0599,   1.9891
c     *,   2.0027,   1.8678,   1.7850,   1.7868,   1.7362
c     *,   1.6478,   1.6254,   1.6166,   1.6478,   1.6516
c     *,   1.5554,   1.5727,   1.5656,   1.5290,   1.5149
c     *,   1.5554,   1.5615,   1.5522,   1.4863,   1.4745
c     *,   1.4745,   1.3991,   1.4139,   1.1686,   1.3458
c     *,   1.3587,   1.2425,   1.3171,   1.2828,   1.2153
c     *,   1.2679,   1.2766,   1.1410,   1.2514,   1.0852
c     *,   1.1445,   1.1672,   1.0624,   1.1089,   1.0899
c     *,   1.0171,   1.0399,   1.0645,   1.0417,    .9934
c     *,   1.0403,   1.0029,   1.0357,   1.0388,   1.0337
c     *,   1.0428,   1.0555,   1.0663,   1.0546,   1.0626
c     *,   1.0013,   1.0705,   1.0711,   1.0852,   1.0830
c     *,   1.1119,   1.0940,   1.0976,   1.0675,   1.1205
c     *,   1.0464,   1.0995,    .9508,   1.0972,   1.1170
c     *,   1.1015,   1.1251,   1.1296,   1.1027,   1.1028
c     *,    .9271,   1.1363,   1.1120,   1.1455,   1.1192
c     *,   1.1498,   1.1491,   1.0108,   1.1282,   1.1550
c     *,   1.1617,   1.1354,   1.1586,   1.1631,   1.1378
c     *,   1.1656,   1.1684,   1.1389,   1.1694,   1.1690
c     *,   1.1702,   1.1705,   1.1390,   1.1714,   1.1716
c     *,   1.1326,   1.1517,   1.1697,   1.1731,   1.1716
c     *,   1.0867,   1.1537,   1.1698,   1.1642,   1.1638
c     *,   1.1199,   1.1635,   1.1713,   1.1631,   1.1601
c     *,   1.1340,   1.0823,   1.1598,   1.1754,   1.1572
c     *,   1.1565,   1.1567,   1.1631,   1.1538,   1.0852
c     *,   1.0342,   1.1645,   1.1452,   1.1099,   1.0940
c     *,   1.1185,   1.1470,   1.1388,   1.1424,   1.0705
c     *,   1.1374,   1.1199,   1.1262,   1.1142,   1.1205
c     *,   1.0867,   1.1095,   1.0734,   1.1234,   1.0925
c     *,   1.1135,   1.1076,   1.1070,   1.0955,   1.1027
c     *,   1.1073,   1.0630,   1.1007,   1.1185,   1.1133
c     *,   1.1128,   1.1028,   1.1033,   1.1086,   1.1120
c     *,   1.1013,   1.0991,   1.1086,   1.1028,   1.1027
c     *,   1.1064,   1.1036,   1.1135,   1.1139,   1.1136
c     *,   1.1145,   1.1166,   1.1152,   1.1168,   1.1172
c     *,   1.1230,   1.1187,   1.1254,   1.1224,   1.1329
c     *,   1.1216
c     */
c      data apnecm/
c     *    2.1288,   2.2007,   2.2500,   2.3044,   2.3529
c     *,   2.9284,   3.5136,   3.6305,   3.7933,   4.1117
c     *,   4.5438,   4.9389,   5.1797,   5.3049,   5.4789
c     *,   5.6476,   6.2773,   6.9855,   7.6283,   8.2211
c     *,   8.7741,   9.2942,   9.7867,  11.5488,  13.7758
c     *,  15.0788,  16.8454,  17.9267,  19.4362,  21.2830
c     *,  22.9819
c     */
c      data apnbmx/
c     *    1.9462,   1.7481,   1.8881,   1.8019,   1.8627
c     *,   1.4745,   1.3231,   1.3762,   1.3351,   1.3505
c     *,   1.2841,   1.3086,   1.2565,   1.3038,   1.2218
c     *,   1.2952,   1.2127,   1.2035,   1.1968,   1.1951
c     *,   1.1604,   1.1794,   1.1747,   1.1724,   1.1595
c     *,   1.1637,   1.1477,   1.1520,   1.1510,   1.1478
c     *,   1.1466
c     */
c      data pipecm/
c     *    1.1050,   1.1154,   1.1165,   1.1256,   1.1273
c     *,   1.1333,   1.1370,   1.1382,   1.1394,   1.1438
c     *,   1.1495,   1.1579,   1.1592,   1.1677,   1.1691
c     *,   1.1697,   1.1757,   1.1771,   1.1777,   1.1784
c     *,   1.1798,   1.1838,   1.1892,   1.1906,   1.1926
c     *,   1.1933,   1.1953,   1.1960,   1.1967,   1.1981
c     *,   1.1994,   1.2008,   1.2015,   1.2022,   1.2028
c     *,   1.2056,   1.2063,   1.2069,   1.2097,   1.2104
c     *,   1.2111,   1.2124,   1.2131,   1.2138,   1.2152
c     *,   1.2166,   1.2172,   1.2193,   1.2200,   1.2214
c     *,   1.2221,   1.2255,   1.2262,   1.2269,   1.2276
c     *,   1.2283,   1.2303,   1.2310,   1.2317,   1.2352
c     *,   1.2358,   1.2365,   1.2372,   1.2400,   1.2407
c     *,   1.2421,   1.2434,   1.2462,   1.2476,   1.2503
c     *,   1.2517,   1.2524,   1.2538,   1.2545,   1.2565
c     *,   1.2613,   1.2627,   1.2696,   1.2716,   1.2730
c     *,   1.2751,   1.2799,   1.2819,   1.2833,   1.2860
c     *,   1.2867,   1.2915,   1.2922,   1.2929,   1.2990
c     *,   1.3010,   1.3024,   1.3119,   1.3126,   1.3180
c     *,   1.3200,   1.3220,   1.3254,   1.3261,   1.3321
c     *,   1.3382,   1.3388,   1.3415,   1.3422,   1.3449
c     *,   1.3469,   1.3522,   1.3608,   1.3615,   1.3622
c     *,   1.3655,   1.3661,   1.3767,   1.3786,   1.3813
c     *,   1.3898,   1.3917,   1.3950,   1.4047,   1.4080
c     *,   1.4112,   1.4157,   1.4163,   1.4176,   1.4208
c     *,   1.4215,   1.4285,   1.4298,   1.4304,   1.4336
c     *,   1.4349,   1.4425,   1.4470,   1.4495,   1.4526
c     *,   1.4539,   1.4652,   1.4677,   1.4702,   1.4764
c     *,   1.4777,   1.4857,   1.4863,   1.4870,   1.4882
c     *,   1.4919,   1.4993,   1.4999,   1.5030,   1.5097
c     *,   1.5103,   1.5121,   1.5146,   1.5212,   1.5249
c     *,   1.5285,   1.5327,   1.5351,   1.5429,   1.5435
c     *,   1.5465,   1.5513,   1.5531,   1.5548,   1.5566
c     *,   1.5637,   1.5655,   1.5673,   1.5708,   1.5720
c     *,   1.5732,   1.5761,   1.5790,   1.5849,   1.5866
c     *,   1.5965,   1.5977,   1.5994,   1.6023,   1.6063
c     *,   1.6121,   1.6133,   1.6144,   1.6264,   1.6287
c     *,   1.6327,   1.6344,   1.6406,   1.6417,   1.6429
c     *,   1.6502,   1.6519,   1.6536,   1.6625,   1.6664
c     *,   1.6670,   1.6676,   1.6687,   1.6692,   1.6714
c     *,   1.6792,   1.6825,   1.6853,   1.6913,   1.6935
c     *,   1.6941,   1.6963,   1.6985,   1.6990,   1.7012
c     *,   1.7170,   1.7175,   1.7208,   1.7240,   1.7267
c     *,   1.7391,   1.7417,   1.7423,   1.7428,   1.7503
c     *,   1.7535,   1.7609,   1.7635,   1.7651,   1.7714
c     *,   1.7735,   1.7767,   1.7793,   1.7798,   1.7908
c     *,   1.7923,   1.7949,   1.7960,   1.8032,   1.8063
c     *,   1.8099,   1.8197,   1.8207,   1.8238,   1.8253
c     *,   1.8289,   1.8304,   1.8320,   1.8386,   1.8436
c     *,   1.8497,   1.8507,   1.8537,   1.8643,   1.8673
c     *,   1.8678,   1.8703,   1.8713,   1.8738,   1.8777
c     *,   1.8792,   1.8812,   1.8827,   1.8881,   1.8911
c     *,   1.8916,   1.8936,   1.8980,   1.9039,   1.9073
c     *,   1.9108,   1.9151,   1.9161,   1.9181,   1.9186
c     *,   1.9249,   1.9283,   1.9317,   1.9331,   1.9350
c     *,   1.9384,   1.9423,   1.9519,   1.9571,   1.9614
c     *,   1.9652,   1.9681,   1.9756,   1.9780,   1.9804
c     *,   1.9893,   1.9921,   1.9968,   1.9987,   2.0034
c     *,   2.0085,   2.0196,   2.0344,   2.0358,   2.0385
c     *,   2.0545,   2.0636,   2.0654,   2.0704,   2.0816
c     *,   2.0839,   2.0933,   2.1106,   2.1151,   2.1257
c     *,   2.1279,   2.1388,   2.1545,   2.1606,   2.1761
c     *,   2.1770,   2.1804,   2.1817,   2.1920,   2.2022
c     *,   2.2060,   2.2187,   2.2480,   2.2646,   2.2934
c     *,   2.3056,   2.3339,   2.3499,   2.3535,   2.3538
c     *,   2.3737,   2.3894,   2.4050,   2.4128,   2.4283
c     *,   2.4398,   2.4513,   2.4703,   2.4892,   2.5080
c     *,   2.5266,   2.5451,   2.5561,   2.5634,   2.5816
c     *,   2.5997,   2.6069,   2.6177,   2.6355,   2.6568
c     *,   2.6708,   2.6987,   2.7057,   2.7195,   2.7401
c     *,   2.7469,   2.7606,   2.7741,   2.7977,   2.8010
c     *,   2.8077,   2.8343,   2.8409,   2.8672,   2.8769
c     *,   2.8997,   2.9093,   2.9350,   2.9414,   2.9636
c     *,   2.9731,   3.0013,   3.0045,   3.0107,   3.0355
c     *,   3.0570,   3.0662,   3.0876,   3.0967,   3.1268
c     *,   3.1328,   3.1566,   3.1862,   3.2067,   3.2155
c     *,   3.2445,   3.2733,   3.3018,   3.3216,   3.3329
c     *,   3.3609,   3.3887,   3.4163,   3.4190,   3.4327
c     *,   3.4436,   3.4707,   3.4869,   3.4976,   3.5244
c     *,   3.5509,   3.6033,   3.6550,   3.7059,   3.7462
c     *,   3.9247,   3.9887,   3.9981,   4.4001,   4.4341
c     *,   4.8193,   4.8387,   5.2120,   5.2246,   5.3889
c     *,   5.5535,   5.5603,   5.7265,   5.8880,   5.8912
c     *,   5.9671,   6.1984,   6.2421,   6.5084,   6.6369
c     *,   6.9138,   7.5617,   8.1584,   8.7144,   8.9794
c     *,   9.2369,   9.7314,   9.9412,  10.2020,  10.6517
c     *,  11.4987,  13.7295,  15.0339,  16.6334,  16.8018
c     *,  17.8835,  18.1439,  21.2400,  22.9386,  24.1342
c     *,  25.2733
c     */
c      data pipbmx/
c     *     .4424,    .5585,    .6180,    .7485,    .7092
c     *,    .8058,    .7777,    .9113,    .8974,    .9934
c     *,   1.0896,   1.2653,   1.3078,   1.5076,   1.5472
c     *,   1.5656,   1.7019,   1.7504,   1.7934,   1.7626
c     *,   1.8168,   1.9706,   2.0163,   2.1140,   2.1749
c     *,   2.0576,   2.1851,   2.1148,   2.1924,   2.1851
c     *,   2.3104,   2.3480,   2.2883,   2.1924,   2.3602
c     *,   2.4570,   2.3262,   2.2708,   2.4722,   2.3296
c     *,   2.5124,   2.3194,   2.4476,   2.4033,   2.5497
c     *,   2.5313,   2.5514,   2.5181,   2.5125,   2.5156
c     *,   2.5105,   2.4398,   2.4201,   2.5074,   2.4978
c     *,   2.4463,   2.4863,   2.4856,   2.4469,   2.5231
c     *,   2.3534,   2.5357,   2.4482,   2.3796,   2.3823
c     *,   2.3582,   2.3650,   2.3870,   2.1705,   2.1185
c     *,   2.2341,   2.2057,   2.2284,   2.1178,   2.1705
c     *,   2.0622,   2.0490,   2.0122,   1.9041,   1.8873
c     *,   1.9091,   1.7897,   1.7499,   1.7535,   1.8797
c     *,   1.8455,   1.6087,   1.6381,   1.6361,   1.6737
c     *,   1.5329,   1.5327,   1.4879,   1.4351,   1.4461
c     *,   1.3854,   1.3517,   1.3334,   1.3695,   1.2741
c     *,   1.2989,   1.2183,   1.2361,   1.2074,   1.1456
c     *,   1.1997,   1.1439,   1.1381,   1.0786,   1.1185
c     *,   1.0693,   1.1339,    .9958,   1.0077,   1.0491
c     *,    .9997,    .9611,    .9452,    .9141,    .9190
c     *,    .8806,    .8885,    .9200,    .8618,    .9680
c     *,    .8702,    .8643,    .8263,    .8234,    .7150
c     *,    .8412,    .8101,    .8704,    .7767,    .8766
c     *,    .7722,    .7694,    .7956,    .7590,    .7128
c     *,    .7159,    .6894,    .6894,    .7139,    .6898
c     *,    .7412,    .6784,    .7181,    .6947,    .7695
c     *,    .6704,    .7110,    .6910,    .6699,    .7199
c     *,    .6935,    .6759,    .6843,    .6933,    .7190
c     *,    .7159,    .6865,    .6958,    .7017,    .6794
c     *,    .7081,    .7004,    .7128,    .7569,    .7269
c     *,    .7356,    .7230,    .7332,    .7548,    .7464
c     *,    .7695,    .7703,    .7707,    .8288,    .8117
c     *,    .8170,    .8054,    .7990,    .7878,    .8556
c     *,    .8423,    .8917,    .8664,    .8253,    .8649
c     *,    .8627,    .8813,    .8461,    .8938,    .9288
c     *,    .8759,    .8345,    .8972,    .8625,    .8818
c     *,    .8983,    .8953,    .9080,    .9277,    .9117
c     *,    .9080,    .8579,    .8974,    .8970,    .9184
c     *,    .9267,    .9184,    .9092,    .8704,    .9146
c     *,    .9184,    .9591,    .8649,    .9296,    .9071
c     *,    .9449,    .9407,    .9536,    .9721,    .9608
c     *,    .9837,    .9233,    .9542,   1.0004,    .9853
c     *,    .9934,   1.0005,   1.0061,    .9358,   1.0180
c     *,   1.0280,   1.0573,    .9982,   1.0357,   1.0555
c     *,    .9877,   1.0682,   1.0694,   1.0817,   1.0720
c     *,   1.1046,   1.1041,   1.0596,   1.1013,   1.1217
c     *,   1.1217,   1.1291,   1.1113,   1.1270,   1.1354
c     *,   1.0799,   1.1199,   1.1439,   1.1459,   1.1213
c     *,   1.0907,   1.1368,   1.1488,   1.1008,   1.1500
c     *,   1.1156,   1.1502,   1.0947,   1.1410,   1.1427
c     *,   1.1480,   1.0729,   1.1634,   1.1319,   1.1480
c     *,   1.1295,   1.1213,   1.0783,   1.1044,   1.1027
c     *,   1.1217,   1.0998,   1.0418,   1.1031,   1.0712
c     *,   1.0720,   1.0600,   1.0540,   1.1180,   1.0490
c     *,   1.0464,   1.0380,   1.0124,   1.0202,    .9821
c     *,   1.0013,    .9945,    .9997,   1.0187,    .9805
c     *,    .9873,    .9796,    .9608,   1.0045,    .9756
c     *,    .9695,    .9674,    .9646,    .9662,    .9659
c     *,    .9654,    .9491,    .9643,    .9566,    .9654
c     *,    .9636,    .9619,    .9676,    .9774,    .9747
c     *,    .9916,    .9843,    .9950,    .9932,    .9641
c     *,    .9895,    .9861,    .9871,    .9905,    .9916
c     *,    .9657,    .9872,    .9770,    .9796,    .9751
c     *,    .9637,    .9689,    .9633,    .9635,    .9601
c     *,    .9580,    .9707,    .9581,    .9540,    .9505
c     *,    .9509,    .9422,    .9494,    .9657,    .9478
c     *,    .9513,    .9541,    .9472,    .9439,    .9640
c     *,    .9452,    .9397,    .9440,    .9390,    .9416
c     *,    .9424,    .9394,    .9334,    .9366,    .9348
c     *,    .9338,    .9414,    .9312,    .9115,    .9286
c     *,    .9373,    .9287,    .9190,    .9248,    .9232
c     *,    .9339,    .9214,    .9201,    .9183,    .9181
c     *,    .9170,    .9150,    .9138,    .9071,    .9117
c     *,    .9106,    .9092,    .9084,    .9253,    .9053
c     *,    .9061,    .9052,    .9132,    .9032,    .9029
c     *,    .9008,    .8973,    .8954,    .8928,    .9115
c     *,    .9021,    .8938,    .8997,    .8907,    .8921
c     *,    .8834,    .8831,    .8795,    .8774,    .8755
c     *,    .8745,    .8682,    .8849,    .8649,    .8705
c     *,    .8616,    .8684,    .8691,    .8634,    .8687
c     *,    .8636,    .8616,    .8581,    .8571,    .8575
c     *,    .8582,    .8551,    .8575,    .8582,    .8618
c     *,    .8597,    .8619,    .8634,    .8693,    .8645
c     *,    .8682,    .8538,    .8759,    .8818,    .8831
c     *,    .8853/
c      data (pimecm(i),i=1,400)/
c     *    1.1046,   1.1133,   1.1394,   1.1425,   1.1495
c     *,   1.1579,   1.1585,   1.1592,   1.1598,   1.1677
c     *,   1.1691,   1.1731,   1.1777,   1.1784,   1.1798
c     *,   1.1831,   1.1858,   1.1879,   1.1892,   1.1906
c     *,   1.1940,   1.1967,   1.1994,   1.2008,   1.2015
c     *,   1.2028,   1.2069,   1.2076,   1.2090,   1.2097
c     *,   1.2111,   1.2124,   1.2131,   1.2159,   1.2166
c     *,   1.2179,   1.2186,   1.2200,   1.2214,   1.2234
c     *,   1.2269,   1.2283,   1.2296,   1.2303,   1.2317
c     *,   1.2324,   1.2352,   1.2358,   1.2365,   1.2372
c     *,   1.2407,   1.2421,   1.2441,   1.2462,   1.2476
c     *,   1.2510,   1.2517,   1.2524,   1.2545,   1.2551
c     *,   1.2579,   1.2586,   1.2593,   1.2607,   1.2613
c     *,   1.2627,   1.2634,   1.2641,   1.2662,   1.2668
c     *,   1.2682,   1.2696,   1.2710,   1.2716,   1.2730
c     *,   1.2758,   1.2778,   1.2806,   1.2819,   1.2826
c     *,   1.2833,   1.2847,   1.2888,   1.2915,   1.2922
c     *,   1.2929,   1.2963,   1.3004,   1.3024,   1.3038
c     *,   1.3058,   1.3065,   1.3072,   1.3112,   1.3119
c     *,   1.3126,   1.3146,   1.3180,   1.3187,   1.3200
c     *,   1.3207,   1.3220,   1.3227,   1.3254,   1.3261
c     *,   1.3301,   1.3308,   1.3321,   1.3335,   1.3362
c     *,   1.3375,   1.3388,   1.3415,   1.3422,   1.3455
c     *,   1.3482,   1.3515,   1.3522,   1.3555,   1.3562
c     *,   1.3615,   1.3628,   1.3641,   1.3655,   1.3681
c     *,   1.3694,   1.3760,   1.3773,   1.3786,   1.3898
c     *,   1.3917,   1.3963,   1.3995,   1.4002,   1.4015
c     *,   1.4047,   1.4099,   1.4112,   1.4163,   1.4176
c     *,   1.4183,   1.4196,   1.4208,   1.4228,   1.4234
c     *,   1.4272,   1.4292,   1.4298,   1.4304,   1.4393
c     *,   1.4425,   1.4470,   1.4476,   1.4489,   1.4507
c     *,   1.4539,   1.4589,   1.4596,   1.4608,   1.4652
c     *,   1.4658,   1.4683,   1.4727,   1.4739,   1.4758
c     *,   1.4764,   1.4777,   1.4795,   1.4808,   1.4833
c     *,   1.4851,   1.4857,   1.4876,   1.4882,   1.4894
c     *,   1.4901,   1.4919,   1.4925,   1.4938,   1.4968
c     *,   1.4974,   1.4981,   1.4993,   1.5030,   1.5054
c     *,   1.5060,   1.5066,   1.5072,   1.5097,   1.5103
c     *,   1.5109,   1.5146,   1.5152,   1.5182,   1.5206
c     *,   1.5212,   1.5224,   1.5231,   1.5249,   1.5309
c     *,   1.5327,   1.5345,   1.5357,   1.5387,   1.5405
c     *,   1.5411,   1.5429,   1.5435,   1.5441,   1.5465
c     *,   1.5489,   1.5507,   1.5513,   1.5519,   1.5548
c     *,   1.5590,   1.5614,   1.5637,   1.5655,   1.5673
c     *,   1.5702,   1.5708,   1.5761,   1.5802,   1.5825
c     *,   1.5831,   1.5843,   1.5849,   1.5866,   1.5890
c     *,   1.5936,   1.5948,   1.5977,   1.5994,   1.6017
c     *,   1.6052,   1.6058,   1.6087,   1.6092,   1.6138
c     *,   1.6173,   1.6178,   1.6190,   1.6236,   1.6258
c     *,   1.6276,   1.6287,   1.6298,   1.6315,   1.6327
c     *,   1.6367,   1.6372,   1.6378,   1.6400,   1.6406
c     *,   1.6423,   1.6480,   1.6513,   1.6519,   1.6525
c     *,   1.6547,   1.6553,   1.6614,   1.6631,   1.6637
c     *,   1.6648,   1.6653,   1.6659,   1.6681,   1.6692
c     *,   1.6703,   1.6720,   1.6726,   1.6731,   1.6742
c     *,   1.6792,   1.6798,   1.6803,   1.6820,   1.6825
c     *,   1.6853,   1.6858,   1.6864,   1.6875,   1.6935
c     *,   1.6957,   1.6963,   1.6979,   1.6990,   1.7067
c     *,   1.7072,   1.7121,   1.7127,   1.7132,   1.7186
c     *,   1.7202,   1.7208,   1.7213,   1.7224,   1.7235
c     *,   1.7240,   1.7267,   1.7278,   1.7294,   1.7342
c     *,   1.7348,   1.7391,   1.7407,   1.7428,   1.7476
c     *,   1.7503,   1.7535,   1.7545,   1.7609,   1.7614
c     *,   1.7667,   1.7714,   1.7735,   1.7741,   1.7767
c     *,   1.7772,   1.7798,   1.7824,   1.7835,   1.7871
c     *,   1.7929,   1.7960,   1.8001,   1.8032,   1.8058
c     *,   1.8130,   1.8182,   1.8212,   1.8218,   1.8238
c     *,   1.8253,   1.8258,   1.8289,   1.8320,   1.8370
c     *,   1.8386,   1.8436,   1.8446,   1.8507,   1.8512
c     *,   1.8557,   1.8588,   1.8638,   1.8678,   1.8713
c     *,   1.8762,   1.8777,   1.8792,   1.8807,   1.8827
c     *,   1.8857,   1.8881,   1.8886,   1.8936,   1.8975
c     *,   1.9010,   1.9063,   1.9132,   1.9186,   1.9205
c     *,   1.9249,   1.9254,   1.9273,   1.9283,   1.9317
c     *,   1.9326,   1.9374,   1.9423,   1.9451,   1.9495
c     *,   1.9519,   1.9538,   1.9614,   1.9709,   1.9733
c     *,   1.9747,   1.9799,   1.9884,   1.9903,   1.9987
c     *,   2.0034,   2.0117,   2.0164,   2.0182,   2.0196
c     *,   2.0330,   2.0335,   2.0358,   2.0540,   2.0581
c     *,   2.0636,   2.0654,   2.0798,   2.0812,   2.0816
c     */
c      data (pimecm(i),i=401,578)/
c     *    2.0933,   2.1013,   2.1075,   2.1221,   2.1261
c     *,   2.1305,   2.1375,   2.1445,   2.1588,   2.1658
c     *,   2.1740,   2.1804,   2.1877,   2.2026,   2.2119
c     *,   2.2166,   2.2187,   2.2305,   2.2510,   2.2626
c     *,   2.2688,   2.2717,   2.3040,   2.3121,   2.3315
c     *,   2.3483,   2.3538,   2.3737,   2.3744,   2.3925
c     *,   2.4050,   2.4105,   2.4128,   2.4267,   2.4302
c     *,   2.4513,   2.4627,   2.4892,   2.5069,   2.5266
c     *,   2.5561,   2.5634,   2.5802,   2.5925,   2.5997
c     *,   2.6284,   2.6355,   2.6557,   2.6708,   2.6987
c     *,   2.7057,   2.7401,   2.7673,   2.7741,   2.7966
c     *,   2.8077,   2.8343,   2.8409,   2.8672,   2.8769
c     *,   2.8997,   2.9093,   2.9318,   2.9341,   2.9414
c     *,   2.9636,   2.9731,   3.0045,   3.0262,   3.0355
c     *,   3.0570,   3.0656,   3.0662,   3.0876,   3.0967
c     *,   3.1268,   3.1566,   3.1774,   3.1862,   3.2067
c     *,   3.2155,   3.2445,   3.2733,   3.3018,   3.3216
c     *,   3.3329,   3.3609,   3.3887,   3.4163,   3.4190
c     *,   3.4436,   3.4675,   3.4707,   3.4869,   3.4976
c     *,   3.5509,   3.6033,   3.6550,   3.6575,   3.7059
c     *,   3.7312,   3.7462,   3.8402,   3.8935,   3.9887
c     *,   4.0004,   4.3873,   4.4341,   4.6811,   4.7210
c     *,   4.7408,   4.8387,   4.8406,   5.0288,   5.0844
c     *,   5.2120,   5.2353,   5.3889,   5.4254,   5.5603
c     *,   5.6123,   5.7265,   5.7787,   5.8880,   5.9451
c     *,   5.9671,   5.9953,   6.0792,   6.1984,   6.2241
c     *,   6.3479,   6.5070,   6.6369,   6.8140,   6.9138
c     *,   7.0734,   7.2450,   7.3962,   7.5617,   7.7092
c     *,   7.9841,   7.9865,   8.1584,   8.1814,   8.3628
c     *,   8.4410,   8.6582,   8.6593,   8.7144,   8.9794
c     *,   9.2369,   9.2845,   9.4874,   9.7314,   9.7775
c     *,   9.9694,  10.2020,  10.3580,  10.4293,  10.5988
c     *,  10.6517,  10.8697,  11.0833,  11.4987,  13.7295
c     *,  15.0339,  16.6334,  16.8018,  17.8835,  18.1439
c     *,  19.3933,  19.6336,  21.2400,  22.9386,  24.1342
c     *,  25.2733,  26.0050,  26.3632
c     */
c      data (pimbmx(i),i=1,400)/
c     *     .5314,    .5314,    .6580,    .7092,    .7506
c     *,    .8215,    .8579,    .8406,    .8349,    .9441
c     *,    .9657,   1.0363,   1.0686,   1.1085,   1.1185
c     *,   1.1899,   1.1658,   1.2218,   1.2322,   1.2653
c     *,   1.3273,   1.2374,   1.3611,   1.3716,   1.3267
c     *,   1.3986,   1.4150,   1.3399,   1.4701,   1.4516
c     *,   1.4723,   1.4955,   1.4494,   1.4161,   1.4558
c     *,   1.4625,   1.5006,   1.5128,   1.5160,   1.4584
c     *,   1.4991,   1.4791,   1.4217,   1.4799,   1.4799
c     *,   1.4588,   1.4844,   1.4172,   1.4273,   1.4439
c     *,   1.4166,   1.4127,   1.3739,   1.3624,   1.3971
c     *,   1.3297,   1.3231,   1.3285,   1.2878,   1.2952
c     *,   1.2679,   1.2641,   1.2976,   1.2386,   1.2463
c     *,   1.3025,   1.2489,   1.2239,   1.1902,   1.2114
c     *,   1.1955,   1.2083,   1.1658,   1.1686,   1.1603
c     *,   1.1424,   1.1185,   1.1160,   1.1013,   1.1070
c     *,   1.1041,   1.0823,   1.0645,   1.0779,   1.0295
c     *,   1.0311,    .9950,   1.0155,    .9950,   1.0029
c     *,    .9767,   1.0249,    .9853,    .9657,    .9906
c     *,    .9608,    .9591,    .9805,    .9458,    .9575
c     *,    .9422,    .9407,    .9558,    .9473,    .9271
c     *,    .9132,    .9236,    .9219,    .9219,    .9575
c     *,    .9097,    .9224,    .9094,    .9202,    .8903
c     *,    .9224,    .9160,    .9122,    .9094,    .8956
c     *,    .9167,    .9271,    .9398,    .9184,    .9277
c     *,    .9310,    .9184,    .9354,    .9219,    .9300
c     *,    .9303,    .9508,    .9451,    .9680,    .9399
c     *,    .9420,    .9387,    .9667,    .9248,    .9915
c     *,    .9698,    .9626,   1.0187,    .9659,    .9441
c     *,    .9766,    .9694,   1.0061,    .9789,    .9757
c     *,    .9803,   1.0495,   1.0218,   1.0026,   1.0706
c     *,   1.0277,   1.0302,   1.0690,   1.0552,   1.0406
c     *,   1.0801,   1.0522,   1.1284,   1.1226,   1.1270
c     *,   1.0934,   1.0968,   1.1381,   1.0915,   1.1543
c     *,   1.1775,   1.2101,   1.1247,   1.1506,   1.1781
c     *,   1.1546,   1.1968,   1.1604,   1.1697,   1.2072
c     *,   1.2058,   1.2489,   1.1912,   1.2141,   1.1928
c     *,   1.1986,   1.1978,   1.2256,   1.2565,   1.2029
c     *,   1.2127,   1.1997,   1.2035,   1.1791,   1.1853
c     *,   1.1914,   1.1927,   1.2616,   1.2303,   1.1848
c     *,   1.1507,   1.2399,   1.1510,   1.1898,   1.1499
c     *,   1.1337,   1.1433,   1.1089,   1.1517,   1.2008
c     *,   1.1598,   1.1030,   1.0798,   1.1165,   1.1213
c     *,   1.1190,   1.0861,   1.0776,   1.0653,   1.0730
c     *,   1.1106,   1.0908,   1.0945,   1.0650,   1.0882
c     *,   1.0720,   1.0540,   1.0925,   1.0718,   1.1377
c     *,   1.0867,   1.0802,   1.1090,   1.0936,   1.0915
c     *,   1.1316,   1.1575,   1.1202,   1.1275,   1.0908
c     *,   1.1142,   1.1507,   1.1890,   1.1863,   1.2114
c     *,   1.1142,   1.2001,   1.1312,   1.1930,   1.2563
c     *,   1.2489,   1.2919,   1.2652,   1.2425,   1.2149
c     *,   1.2386,   1.2616,   1.3207,   1.3299,   1.2816
c     *,   1.3296,   1.3562,   1.2957,   1.3831,   1.3028
c     *,   1.3806,   1.3704,   1.3202,   1.4082,   1.3482
c     *,   1.3730,   1.3351,   1.3886,   1.3348,   1.3957
c     *,   1.3929,   1.3790,   1.3791,   1.3946,   1.3641
c     *,   1.3877,   1.3566,   1.3564,   1.3762,   1.3253
c     *,   1.3739,   1.3252,   1.3234,   1.3423,   1.3159
c     *,   1.2643,   1.2700,   1.2584,   1.2466,   1.2841
c     *,   1.2191,   1.2052,   1.2794,   1.2835,   1.2374
c     *,   1.2061,   1.1936,   1.1947,   1.2476,   1.1481
c     *,   1.1496,   1.2187,   1.1427,   1.2101,   1.1192
c     *,   1.1445,   1.1371,   1.1613,   1.0994,   1.1020
c     *,   1.1233,   1.1041,   1.0896,   1.0890,   1.0880
c     *,   1.0802,   1.0946,   1.0940,   1.0670,   1.0761
c     *,   1.0817,   1.0702,   1.0723,   1.0779,   1.0789
c     *,   1.0733,   1.0771,   1.0767,   1.0710,   1.0705
c     *,   1.0720,   1.0729,   1.0472,   1.0959,   1.0788
c     *,   1.0739,   1.0645,   1.0681,   1.0749,   1.0789
c     *,   1.0799,   1.0795,   1.0763,   1.0588,   1.0805
c     *,   1.0724,   1.0660,   1.0663,   1.0801,   1.0988
c     *,   1.0797,   1.0816,   1.0687,   1.0690,   1.0724
c     *,   1.0721,   1.0777,   1.0650,   1.0715,   1.0617
c     *,   1.0690,   1.0678,   1.0714,   1.0323,   1.0672
c     *,   1.0630,   1.0564,   1.0600,   1.0636,   1.0499
c     *,   1.0479,   1.0481,   1.0463,   1.0510,   1.0450
c     *,   1.0171,   1.0517,   1.0498,   1.0510,   1.0412
c     *,   1.0374,   1.0501,   1.0388,   1.0449,   1.0449
c     *,   1.0311,   1.0510,   1.0505,   1.0449,   1.0567
c     *,   1.0572,   1.0495,   1.0632,   1.0461,   1.0403
c     */
c      data (pimbmx(i),i=401,578)/
c     *    1.0629,   1.0642,   1.0525,   1.0709,   1.0585
c     *,   1.0612,   1.0717,   1.0731,   1.0660,   1.0761
c     *,   1.0696,   1.0767,   1.0755,   1.0720,   1.0615
c     *,   1.0665,   1.0660,   1.0714,   1.0670,   1.0650
c     *,   1.0627,   1.0621,   1.0499,   1.0499,   1.0457
c     *,   1.0405,   1.0372,   1.0412,   1.0331,   1.0299
c     *,   1.0299,   1.0241,   1.0321,   1.0306,   1.0232
c     *,   1.0252,   1.0171,   1.0211,   1.0226,   1.0178
c     *,   1.0083,   1.0157,   1.0138,    .9918,   1.0140
c     *,   1.0028,   1.0120,   1.0077,   1.0088,    .9979
c     *,   1.0052,   1.0025,    .9876,    .9987,    .9990
c     *,    .9944,    .9845,    .9918,    .9777,    .9892
c     *,    .9754,    .9856,    .9812,    .9816,    .9831
c     *,    .9674,    .9800,    .9781,    .9679,    .9756
c     *,    .9805,    .9759,    .9739,    .9672,    .9718
c     *,    .9701,    .9682,    .9632,    .9664,    .9538
c     *,    .9647,    .9628,    .9613,    .9596,    .9449
c     *,    .9585,    .9572,    .9553,    .9540,    .9624
c     *,    .9527,    .9468,    .9513,    .9525,    .9503
c     *,    .9479,    .9459,    .9438,    .9424,    .9422
c     *,    .9407,    .9457,    .9399,    .9385,    .9368
c     *,    .9409,    .9248,    .9209,    .9034,    .8974
c     *,    .9150,    .9089,    .9145,    .9146,    .9080
c     *,    .9044,    .9082,    .9069,    .9062,    .8938
c     *,    .9034,    .9045,    .9011,    .8921,    .8979
c     *,    .8925,    .8982,    .8975,    .8976,    .8947
c     *,    .8962,    .8932,    .8916,    .8913,    .8888
c     *,    .8889,    .8892,    .8880,    .8881,    .8842
c     *,    .8815,    .8826,    .8832,    .8813,    .8826
c     *,    .8791,    .8835,    .8831,    .8833,    .8782
c     *,    .8789,    .8813,    .8797,    .8780,    .8744
c     *,    .8782,    .8842,    .8777,    .8777,    .8789
c     *,    .8815,    .8800,    .8839,    .8740,    .8732
c     *,    .8751,    .8784,    .8757,    .8779,    .8630
c     *,    .8802,    .8775,    .8851,    .8874,    .8903
c     *,    .8935,    .8965,    .8965
c     */
c      data kmpecm/
c     *    1.4691,   1.4720,   1.4750,   1.4780,   1.4811
c     *,   1.4837,   1.4843,   1.4860,   1.4876,   1.4910
c     *,   1.4944,   1.4979,   1.5014,   1.5032,   1.5050
c     *,   1.5087,   1.5091,   1.5124,   1.5132,   1.5162
c     *,   1.5170,   1.5200,   1.5220,   1.5239,   1.5278
c     *,   1.5318,   1.5354,   1.5358,   1.5362,   1.5378
c     *,   1.5399,   1.5440,   1.5523,   1.5607,   1.5654
c     *,   1.5688,   1.5775,   1.5784,   1.5863,   1.5916
c     *,   1.5947,   1.6023,   1.6050,   1.6055,   1.6086
c     *,   1.6145,   1.6159,   1.6172,   1.6191,   1.6236
c     *,   1.6319,   1.6328,   1.6332,   1.6420,   1.6461
c     *,   1.6466,   1.6522,   1.6563,   1.6582,   1.6614
c     *,   1.6642,   1.6694,   1.6712,   1.6717,   1.6768
c     *,   1.6806,   1.6811,   1.6839,   1.6843,   1.6867
c     *,   1.6885,   1.6961,   1.6965,   1.7003,   1.7022
c     *,   1.7083,   1.7087,   1.7172,   1.7177,   1.7181
c     *,   1.7229,   1.7243,   1.7276,   1.7342,   1.7374
c     *,   1.7436,   1.7459,   1.7483,   1.7539,   1.7610
c     *,   1.7629,   1.7633,   1.7718,   1.7770,   1.7789
c     *,   1.7793,   1.7817,   1.7840,   1.7873,   1.7892
c     *,   1.8028,   1.8037,   1.8075,   1.8135,   1.8140
c     *,   1.8219,   1.8247,   1.8261,   1.8308,   1.8369
c     *,   1.8401,   1.8406,   1.8410,   1.8475,   1.8480
c     *,   1.8489,   1.8540,   1.8559,   1.8605,   1.8647
c     *,   1.8721,   1.8744,   1.8767,   1.8785,   1.8836
c     *,   1.8951,   1.8955,   1.8983,   1.9001,   1.9029
c     *,   1.9065,   1.9184,   1.9202,   1.9243,   1.9348
c     *,   1.9393,   1.9411,   1.9434,   1.9483,   1.9547
c     *,   1.9628,   1.9637,   1.9659,   1.9699,   1.9798
c     *,   1.9838,   1.9923,   1.9958,   2.0047,   2.0127
c     *,   2.0149,   2.0162,   2.0255,   2.0277,   2.0308
c     *,   2.0417,   2.0430,   2.0487,   2.0579,   2.0653
c     *,   2.0679,   2.0713,   2.0813,   2.0882,   2.0925
c     *,   2.1028,   2.1105,   2.1126,   2.1211,   2.1232
c     *,   2.1258,   2.1351,   2.1444,   2.1507,   2.1528
c     *,   2.1654,   2.1675,   2.1687,   2.1737,   2.1837
c     *,   2.1862,   2.1937,   2.2044,   2.2131,   2.2143
c     *,   2.2274,   2.2356,   2.2478,   2.2546,   2.2659
c     *,   2.2756,   2.2836,   2.2976,   2.2996,   2.3162
c     *,   2.3166,   2.3296,   2.3335,   2.3363,   2.3535
c     *,   2.3562,   2.3725,   2.3729,   2.3748,   2.3887
c     *,   2.3918,   2.3933,   2.4006,   2.4109,   2.4174
c     *,   2.4223,   2.4299,   2.4352,   2.4488,   2.4518
c     *,   2.4675,   2.4705,   2.4861,   2.4887,   2.4936
c     *,   2.5046,   2.5230,   2.5412,   2.5593,   2.5701
c     *,   2.5773,   2.5952,   2.6023,   2.6059,   2.6130
c     *,   2.6306,   2.6447,   2.6482,   2.6656,   2.6795
c     *,   2.6829,   2.7002,   2.7173,   2.7848,   2.8540
c     *,   2.9216,   2.9407,   2.9878,   3.0096,   3.0250
c     *,   3.0526,   3.1783,   3.2191,   3.2993,   3.3887
c     *,   3.5239,   3.6924,   3.7800,   3.9967,   4.0200
c     *,   4.1348,   4.4617,   4.4826,   4.7663,   4.8636
c     *,   4.9779,   5.1080,   5.1263,   5.2349,   5.3237
c     *,   5.4110,   5.5479,   5.5816,   5.8281,   5.9080
c     *,   6.0801,   6.2173,   6.3664,   6.6545,   6.9306
c     *,   7.2610,   7.5770,   7.7605,   7.8207,   7.9985
c     *,   8.1725,   8.2297,   8.4546,   8.7275,   8.9922
c     *,   9.2493,   9.4994,   9.7431,   9.9809,  10.2131
c     *,  11.5086,  13.7378,  15.0415,  16.8085,  17.8898
c     *,  19.3991,  21.2453,  22.9435,  24.1389
c     */
c      data kmpbmx/
c     *    1.9033,   1.7662,   1.7298,   1.7544,   1.5461
c     *,   1.6991,   1.6205,   1.5898,   1.5817,   1.5023
c     *,   1.5554,   1.5086,   1.5065,   1.4948,   1.4799
c     *,   1.4927,   1.4854,   1.6136,   1.6017,   1.7312
c     *,   1.4884,   1.7075,   1.5574,   1.5260,   1.5002
c     *,   1.4571,   1.3991,   1.3219,   1.3434,   1.3635
c     *,   1.3315,   1.2966,   1.2213,   1.1604,   1.1930
c     *,   1.1382,   1.1354,   1.1583,   1.1185,   1.1156
c     *,   1.0645,   1.0247,   1.0969,   1.0779,   1.1083
c     *,   1.0749,   1.0155,   1.0464,   1.0754,   1.0202
c     *,   1.0470,   1.0615,   1.0540,   1.0357,   1.0605
c     *,   1.0428,   1.0265,   1.0187,   1.0611,   1.0615
c     *,   1.0295,   1.0651,   1.0823,   1.0077,   1.0764
c     *,   1.1298,   1.1013,   1.1080,   1.1070,   1.1368
c     *,   1.1382,   1.1410,   1.1418,   1.1312,   1.1466
c     *,   1.1378,   1.1241,   1.0896,   1.1264,   1.0720
c     *,   1.1368,   1.1013,   1.1442,   1.1466,   1.1612
c     *,   1.1726,   1.1755,   1.1576,   1.1726,   1.1848
c     *,   1.1617,   1.2029,   1.1686,   1.2274,   1.2252
c     *,   1.1726,   1.2256,   1.2244,   1.2008,   1.2332
c     *,   1.2828,   1.2548,   1.2412,   1.2887,   1.2853
c     *,   1.2653,   1.2540,   1.2527,   1.2449,   1.2118
c     *,   1.1781,   1.1902,   1.1767,   1.1859,   1.1410
c     *,   1.1594,   1.1800,   1.1185,   1.1256,   1.1354
c     *,   1.1131,   1.1562,   1.1143,   1.1382,   1.0841
c     *,   1.0587,   1.0632,   1.0311,   1.0617,   1.0150
c     *,   1.0309,   1.0174,   1.0110,   1.0171,    .9961
c     *,    .9719,    .9938,    .9869,    .9953,    .9966
c     *,   1.0189,    .9977,    .9918,    .9948,   1.0034
c     *,    .9737,   1.0066,   1.0137,   1.0041,   1.0171
c     *,   1.0303,   1.0223,   1.0322,    .9751,   1.0331
c     *,   1.0118,   1.0403,   1.0399,   1.0429,   1.0171
c     *,   1.0279,   1.0467,   1.0414,   1.0217,   1.0432
c     *,   1.0379,   1.0280,   1.0351,   1.0171,   1.0278
c     *,   1.0314,   1.0213,   1.0133,   1.0240,   1.0080
c     *,   1.0018,    .9964,    .9984,    .9889,    .9910
c     *,    .9903,    .9801,    .9852,    .9853,    .9847
c     *,    .9800,    .9832,    .9770,    .9749,    .9754
c     *,    .9723,    .9741,    .9744,    .9738,    .9751
c     *,    .9780,    .9738,    .9738,    .9684,    .9712
c     *,    .9669,    .9680,    .9671,    .9576,    .9619
c     *,    .9624,    .9585,    .9588,    .9586,    .9525
c     *,    .9253,    .9518,    .9503,    .9491,    .9463
c     *,    .9476,    .9420,    .9458,    .9434,    .9341
c     *,    .9444,    .9412,    .9393,    .9395,    .9210
c     *,    .9370,    .9358,    .8974,    .9400,    .9342
c     *,    .9305,    .9141,    .9271,    .9267,    .9228
c     *,    .9233,    .9219,    .9247,    .9296,    .9253
c     *,    .9089,    .8992,    .8946,    .8992,    .9062
c     *,    .9069,    .8874,    .9219,    .8746,    .8795
c     *,    .8740,    .8704,    .8785,    .8795,    .8667
c     *,    .8849,    .8510,    .8463,    .8612,    .8292
c     *,    .8292,    .8387,    .8273,    .8273,    .8292
c     *,    .8292,    .8349,    .8234,    .8349,    .8176
c     *,    .8292,    .8180,    .8193,    .8154,    .8130
c     *,    .8121,    .8145,    .8078,    .7959,    .8088
c     *,    .8075,    .8064,    .8056,    .8086,    .8048
c     *,    .8080,    .8068,    .8050,    .8042,    .8050
c     *,    .8054,    .8064,    .8096,    .8095,    .8107
c     *,    .8135,    .8234,    .8238,    .8263
c     */
c      data kmnecm/
c     *    1.6212,   1.6781,   1.7053,   1.7863,   1.7882
c     *,   1.8271,   1.9025,   2.0152,   2.1237,   2.2157
c     *,   2.4252,   2.6054,   2.6160,   2.9441,   3.5279
c     *,   3.6965,   4.0245,   4.4666,   4.8690,   5.1136
c     *,   5.2406,   5.4169,   5.5877,   5.9144,   6.2241
c     *,   6.9381,   7.5852,   8.1813,   8.7369,   9.2592
c     *,   9.7536,  10.2241,  11.5209,  13.7525,  15.0575
c     *,  16.8264,  17.9089,  19.4198,  21.2680,  22.9680
c     *,  24.1646
c     */
c      data kmnbmx/
c     *     .9132,    .9966,    .9772,    .9966,   1.1382
c     *,   1.0794,    .9674,    .9167,    .8795,    .8500
c     *,    .8482,    .8444,    .8330,    .8078,    .8349
c     *,    .8368,    .7919,    .8146,    .8019,    .8156
c     *,    .7999,    .7961,    .8038,    .8038,    .7949
c     *,    .7868,    .7864,    .7910,    .7850,    .7870
c     *,    .7906,    .7885,    .7945,    .7971,    .8011
c     *,    .8003,    .8029,    .8082,    .8088,    .8156
c     *,    .8189
c     */
c      data kppecm/
c     *    1.4447,   1.4516,   1.4522,   1.4585,   1.4663
c     *,   1.4750,   1.5036,   1.5050,   1.5091,   1.5239
c     *,   1.5378,   1.5523,   1.5654,   1.5784,   1.5916
c     *,   1.5929,   1.6014,   1.6037,   1.6050,   1.6150
c     *,   1.6159,   1.6191,   1.6259,   1.6264,   1.6268
c     *,   1.6328,   1.6378,   1.6461,   1.6517,   1.6587
c     *,   1.6605,   1.6652,   1.6792,   1.6843,   1.6853
c     *,   1.6928,   1.7073,   1.7101,   1.7210,   1.7294
c     *,   1.7374,   1.7422,   1.7483,   1.7539,   1.7643
c     *,   1.7662,   1.7704,   1.7789,   1.7793,   1.7864
c     *,   1.7897,   1.8028,   1.8070,   1.8135,   1.8191
c     *,   1.8327,   1.8355,   1.8373,   1.8517,   1.8587
c     *,   1.8605,   1.8679,   1.8725,   1.8813,   1.8836
c     *,   1.9038,   1.9070,   1.9289,   1.9298,   1.9321
c     *,   1.9524,   1.9533,   1.9749,   1.9807,   1.9950
c     *,   1.9972,   1.9994,   2.0074,   2.0193,   2.0435
c     *,   2.0492,   2.0635,   2.0653,   2.0851,   2.1040
c     *,   2.1066,   2.1083,   2.1279,   2.1296,   2.1490
c     *,   2.1507,   2.1717,   2.1908,   2.1924,   2.2110
c     *,   2.2131,   2.2213,   2.2319,   2.2335,   2.2538
c     *,   2.2724,   2.2740,   2.2940,   2.3123,   2.3138
c     *,   2.3375,   2.3531,   2.3725,   2.3903,   2.3918
c     *,   2.4109,   2.4197,   2.4299,   2.4413,   2.4488
c     *,   2.4675,   2.4861,   2.5046,   2.5230,   2.5266
c     *,   2.5412,   2.5521,   2.5593,   2.5773,   2.5952
c     *,   2.6130,   2.6306,   2.6482,   2.6656,   2.6829
c     *,   2.7002,   2.7173,   3.0096,   3.0556,   3.1754
c     *,   3.3887,   3.5239,   3.7800,   4.0200,   4.1348
c     *,   4.4617,   4.6469,   4.7663,   4.8636,   4.9591
c     *,   5.1263,   5.2349,   5.4110,   5.5816,   5.7308
c     *,   5.9080,   6.0646,   6.2173,   6.9306,   7.5770
c     *,   7.8207,   8.1725,   8.7275,   9.2493,   9.7431
c     *,  10.2131,  11.5086,  13.7378,  15.0415,  16.6402
c     *,  16.8085,  17.8898,  18.1501,  19.3991,  21.2453
c     *,  22.9435,  24.1389
c     */
c      data kppbmx/
c     *     .5412,    .6308,    .6024,    .6050,    .5971
c     *,    .6037,    .6232,    .6103,    .6482,    .6601
c     *,    .6386,    .6466,    .6438,    .6204,    .6482
c     *,    .6358,    .6333,    .6445,    .6443,    .6346
c     *,    .6410,    .6227,    .6283,    .6308,    .6403
c     *,    .6290,    .6457,    .5984,    .6333,    .5955
c     *,    .5944,    .6295,    .6346,    .6090,    .6433
c     *,    .6383,    .6482,    .6425,    .6543,    .6588
c     *,    .6652,    .6768,    .6730,    .6723,    .6815
c     *,    .7040,    .6898,    .7014,    .7001,    .7181
c     *,    .7130,    .7159,    .7067,    .7440,    .7345
c     *,    .7365,    .7485,    .7382,    .7474,    .7574
c     *,    .7588,    .7559,    .7590,    .7582,    .7668
c     *,    .7592,    .7682,    .7661,    .7697,    .7548
c     *,    .7661,    .7626,    .7626,    .7563,    .7590
c     *,    .7578,    .7611,    .7557,    .7555,    .7506
c     *,    .7498,    .7517,    .7508,    .7540,    .7464
c     *,    .7538,    .7512,    .7527,    .7534,    .7527
c     *,    .7565,    .7521,    .7529,    .7525,    .7444
c     *,    .7517,    .7334,    .7485,    .7491,    .7510
c     *,    .7466,    .7476,    .7478,    .7472,    .7485
c     *,    .7378,    .7451,    .7468,    .7474,    .7476
c     *,    .7459,    .7410,    .7461,    .7457,    .7414
c     *,    .7464,    .7457,    .7444,    .7444,    .7444
c     *,    .7442,    .7291,    .7421,    .7429,    .7421
c     *,    .7397,    .7386,    .7373,    .7389,    .7384
c     *,    .7384,    .7386,    .7378,    .7569,    .7878
c     *,    .7548,    .7356,    .7526,    .7421,    .7715
c     *,    .7568,    .7590,    .7777,    .7421,    .7632
c     *,    .7464,    .7442,    .7548,    .7367,    .7736
c     *,    .7378,    .7421,    .7455,    .7502,    .7510
c     *,    .7653,    .7529,    .7580,    .7544,    .7614
c     *,    .7605,    .7661,    .7703,    .7805,    .7883
c     *,    .7850,    .7907,    .7611,    .7961,    .8023
c     *,    .8068,    .8111
c     */
c      data kpnecm/
c     *    1.6875,   1.7430,   1.7670,   1.7816,   1.7905
c     *,   1.8144,   1.8383,   1.8615,   1.8749,   1.8846
c     *,   1.9080,   1.9308,   1.9345,   1.9535,   1.9760
c     *,   1.9974,   1.9983,   2.0205,   2.0460,   2.0647
c     *,   2.0678,   2.0864,   2.1066,   2.1079,   2.1109
c     *,   2.1292,   2.1322,   2.1504,   2.1533,   2.1743
c     *,   2.1922,   2.1951,   2.2157,   2.2239,   2.2334
c     *,   2.2362,   2.2565,   2.2739,   2.2767,   2.2967
c     *,   2.3138,   2.3166,   2.3402,   2.3559,   2.3753
c     *,   2.3919,   2.3946,   2.4138,   2.4328,   2.4517
c     *,   2.4704,   2.4891,   2.5076,   2.5259,   2.5442
c     *,   2.5551,   2.5623,   2.5803,   2.5982,   2.6160
c     *,   2.6337,   2.6513,   2.6687,   2.6861,   2.7033
c     *,   2.7204,   3.5279,   4.0245,   4.4666,   4.8690
c     *,   5.2406,   5.4169,   5.5877,   5.9144,   6.2241
c     *,   6.9381,   7.5852,   8.1813,   8.7369,   9.2592
c     *,   9.7536,  10.2241,  11.5209,  13.7525,  15.0575
c     *,  16.8264,  17.9089,  19.4198,  21.2680,  22.9680
c     *,  24.1646
c     */
c      data kpnbmx/
c     *     .7024,    .7324,    .7485,    .7527,    .7680
c     *,    .7758,    .8100,    .8224,    .7611,    .8151
c     *,    .8031,    .7915,    .7674,    .7842,    .7822
c     *,    .7590,    .7791,    .7767,    .7758,    .7734
c     *,    .7754,    .7709,    .7674,    .7713,    .7742
c     *,    .7752,    .7748,    .7721,    .7680,    .7707
c     *,    .7674,    .7713,    .7715,    .7695,    .7684
c     *,    .7734,    .7682,    .7709,    .7672,    .7659
c     *,    .7653,    .7653,    .7506,    .7626,    .7624
c     *,    .7701,    .7588,    .7622,    .7592,    .7491
c     *,    .7588,    .7574,    .7592,    .7582,    .7571
c     *,    .7464,    .7559,    .7538,    .7529,    .7529
c     *,    .7534,    .7538,    .7487,    .7487,    .7498
c     *,    .7474,    .7464,    .7485,    .7464,    .7485
c     *,    .7464,    .7565,    .7442,    .7485,    .7545
c     *,    .7555,    .7542,    .7634,    .7653,    .7686
c     *,    .7671,    .7723,    .7695,    .7785,    .7824
c     *,    .7905,    .7927,    .7923,    .8052,    .8100
c     *,    .8137
c     */
cc      data lpecm/
cc     *    2.0577,   2.0583,   2.0593,   2.0595,   2.0609
cc     *,   2.0617,   2.0629,   2.0642,   2.0643,   2.0647
cc     *,   2.0666,   2.0671,   2.0683,   2.0692,   2.0709
cc     *,   2.0720,   2.0783,   2.0818,   2.0935,   2.1022
cc     *,   2.1058,   2.1117,   2.1326,   2.1559,   2.1731
cc     *,   2.1811,   2.2079,   2.2360,   2.2505,   2.3107
cc     *,   2.3733,   2.4534,   2.4695,   2.9278,   5.2476
cc     */
cc      data lpbmx/
cc     *    2.5793,   2.3937,   1.9381,   2.3736,   2.0342
cc     *,   2.2068,   1.9381,   1.8797,   1.7841,   1.7930
cc     *,   1.2676,   1.6641,   1.4604,   1.0403,   1.3470
cc     *,   1.0420,    .8722,    .8368,    .5323,    .5352
cc     *,    .8867,    .5382,    .8106,    .8043,    .8058
cc     *,    .8630,    .7031,    .7777,    .9441,   1.0403
cc     *,    .9772,    .6628,   1.1968,   1.0555,   1.0495
cc     */
c      data pi1ecm/
c     *     .45000,    .55000,    .62500,    .68750,    .73750
c     *,    .78750,    .83750,    .88750,    .95000,   1.02500
c     *,   1.10000,   1.17500
c     */
c      data pi1bmx/
c     *     .48533,    .47873,    .55852,    .58087,    .55279
c     *,    .49185,    .47203,    .46181,    .46181,    .47203
c     *,    .41459,    .40684
c     */
c      data pi2ecm/
c     *     .45000,    .55000,    .62500,    .68750,    .73750
c     *,    .78750,    .83750,    .88750,    .95000,   1.02500
c     *,   1.10000,   1.17500
c     */
c      data pi2bmx/
c     *     .45135,    .80385,    .87767,   1.28779,   1.81771
c     *,   2.07296,   1.55536,    .95912,    .94407,    .77768
c     *,    .70467,    .54408
c     */
c      data pi3ecm/
c     *     .45000,    .55000,    .61250,    .65000,    .67500
c     *,    .70000,    .72500,    .75000,    .77500,    .80000
c     *,    .82500,    .85000,    .87500,    .90000,    .93750
c     *,   1.02500,   1.10000,   1.17500
c     */
c      data pi3bmx/
c     *     .50463,    .79988,   1.13541,   1.30497,   1.34462
c     *,   1.52957,   1.63809,   1.94216,   1.98511,   1.66316
c     *,   1.69163,   1.62933,   1.06152,    .96574,    .83302
c     *,    .71365,    .66517,    .59173
c     */
c      data pi4ecm/
c     *     .28500,    .29500,    .30500,    .31500,    .37000
c     *,    .46000,    .54375,    .61250,    .65000,    .67500
c     *,    .70000,    .72500,    .75000,    .77500,    .80000
c     *,    .82500,    .85000,    .87500,    .90125,    .92750
c     *,    .96000
c     */
c      data pi4bmx/
c     *     .29854,    .46181,    .54115,    .54115,    .51709
c     *,    .72031,    .80976,    .96078,    .90798,   1.07788
c     *,   1.37389,   1.70101,   2.02716,   1.83687,   1.52644
c     *,   1.57469,   1.35640,    .93390,    .98531,    .84062
c     *,    .79188
c     */
c      data pi5ecm/0.75,0.80,0.85,0.90,0.95,1.00,1.05,1.10,1.15/
c      data pi5bmx/0.4,0.5,0.7,0.9,1.1,1.1,0.7,0.6,0.5/
c
c      call utpri('jintcs',ish,ishini,7)
c      
c      ics=0
c
c      if(idptl(i).gt.10000)then
c       call idquad(i,nqi,nai,jci)
c       if(nqi.eq.0)then
c        idi=999
c       else
c        idi=9999
c       endif
c      else
c       idi=idptl(i)
c      endif
c
c      if(idptl(j).gt.10000)then
c       call idquad(j,nqj,naj,jcj)
c       if(nqj.eq.0)then
c        idj=999
c       else
c        idj=9999
c       endif
c      else
c       idj=idptl(j)
c      endif
c
c      do k=1,nflav
c       kc(k)=jc(k,1)-jc(k,2)
c       if(nq.lt.0)kc(k)=-kc(k)
c      enddo
c      if(nq.lt.0)nq=-nq
c
c      if(ish.ge.7)then
c       write(ifch,*)'i:',i,' id_i:',idptl(i),' j:',j,' id_j:',idptl(j)
c       write(ifch,*)'E_cm:',ecm,' jc:',(jc(k,1),k=1,6),(jc(k,2),k=1,6)
c       write(ifch,*)'nq:',nq,' kc:',(kc(k),k=1,6)
c      endif
c
cc check minimal kinetic energy
c      ekin=ecm-pptl(5,i)-pptl(5,j)
c      if(ekin.lt.amimel)goto1000
c
cc -------------------------------------------------------------------------
c      if(iabs(idi).gt.1000.or.iabs(idj).gt.1000)then   !baryon involved
cc -------------------------------------------------------------------------
c      
c       if(nq.eq.6)then !------------baryon-baryon -----> 
c          
c        if(kc(1).eq.kc(2))then    !isospin_z=0
cc pn
c         if(ish.ge.7)write(ifch,*)'sig_pn chosen'
c         if(ecm.lt.2)then
c          call utindx(npn,pnecm,ecm,iecm)
c          bmx=pnbmx(iecm)
c         else
c          p=ecm*sqrt((ecm**2/4./.94**2)-1.)
c          if(p.lt.1.)then
c           sig=33.+196.*(abs(0.95-p))**2.5
c          elseif(p.lt.2)then
c           sig=24.2+8.9*p
c          else
c           sig=47.3+.513*(alog(p)**2)-4.27*alog(p)
c          endif
c          bmx=sqrt(sig/10./pi)
c         endif
c        else                      !isospin_z=+-1
cc pp
c         if(ish.ge.7)write(ifch,*)'sig_pp chosen'
c         if(ecm.le.1.882)then
c          bmx=ppbmx(1) 
c          if(ish.ge.7)write(ifch,*)'b_mx:',bmx
c         elseif(ecm.le.1.887)then
c          bmx=ppbmx(2)
c          if(ish.ge.7)write(ifch,*)'b_mx:',bmx
c         elseif(ecm.le.1.893)then
c          bmx=ppbmx(3)
c          if(ish.ge.7)write(ifch,*)'b_mx:',bmx
c         elseif(ecm.le.1.899)then
c          bmx=ppbmx(4)
c         elseif(ecm.le.1.909)then
c          bmx=ppbmx(5)
c         elseif(ecm.le.1.913)then
c          bmx=ppbmx(6)
c         elseif(ecm.le.1.918)then
c          bmx=ppbmx(7)
c         else
c          p=ecm*sqrt((ecm**2/4./.94**2)-1.)
c          if(p.lt.0.8)then
c           sig=23.5+1000.*(0.7-p)**4
c          elseif(p.lt.1.5)then
c           sig=23.5+24.6/(1+exp(-(p-1.2)/0.1))
c          elseif(p.lt.5)then
c           sig=41.+60.*(p-0.9)*exp(-1.2*p)
c          else
c           sig=48+.522*(alog(p)**2)-4.51*alog(p)
c          endif
c          bmx=sqrt(sig/10./pi)
c         endif
c        endif
c        
c       elseif(nq.eq.0)then !-----------baryon-antibaryon ---->  
c          
c        if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0)then  !p-ap, n-an
cc app
c         if(ish.ge.7)write(ifch,*)'sig_app chosen'
c         if(ecm.lt.2)then
c          call utindx(napp,appecm,ecm,iecm)
c          bmx=appbmx(iecm)
c         else
c          p=ecm*sqrt((ecm**2/4./.94**2)-1.)
c          sig=38.4+77.6*p**(-.64)+.26*(alog(p)**2)-1.2*alog(p)
c          bmx=sqrt(sig/10./pi)
c         endif
c        else                                              !p-an, n-ap
cc apn
c         if(ish.ge.7)write(ifch,*)'sig_apn chosen'
c         if(ecm.lt.2)then
c          call utindx(napn,apnecm,ecm,iecm)
c          bmx=apnbmx(iecm)
c         else
c          p=ecm*sqrt((ecm**2/4./.94**2)-1.)
c          sig=133.6*p**(-.7)-1.22*(alog(p)**2)+13.7*alog(p)
c          bmx=sqrt(sig/10./pi)
c         endif
c        endif
c        
c       elseif(nq.eq.3)then !----------baryon-meson ----> 
c                          
c        if(kc(3).eq.0)then        
cc no kaons involved(except for K-L)
c         if(kc(1).eq.0.or.kc(2).eq.0)then
cc pip
c          if(ish.ge.7)write(ifch,*)'sig_pip chosen'
c          if(ecm.lt.2.5)then
c           call utindx(npip,pipecm,ecm,iecm)
c           bmx=pipbmx(iecm)
c          else
c           p=sqrt(((ecm**2-0.94**2-0.14**2)/2./0.94)**2.-0.14**2)
c           sig=16.4+19.3*p**(-.42)+.19*(alog(p)**2)
c           bmx=sqrt(sig/10./pi)
c          endif
c         else
cc pim
c          if(ish.ge.7)write(ifch,*)'sig_pim chosen'
c          if(ecm.lt.2.5)then
c           call utindx(npim,pimecm,ecm,iecm)
c           bmx=pimbmx(iecm)
c          else
c           p=sqrt(((ecm**2-0.94**2-0.14**2)/2./0.94)**2.-0.14**2)
c           sig=33.0+14.0*p**(-1.36)+.456*(alog(p)**2)-4.03*alog(p)
c           bmx=sqrt(sig/10./pi)
c          endif
c         endif
c        elseif(kc(3).eq.1)then      
cc strange particles involved
c         if(kc(1).eq.0.or.kc(2).eq.0)then
cc kmn
c          if(ish.ge.7)write(ifch,*)'sig_kmn chosen'
c          if(ecm.lt.2.5)then
c           call utindx(nkmn,kmnecm,ecm,iecm)
c           bmx=kmnbmx(iecm)
c          else
c           p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2)
c           sig=25.2+.38*(alog(p)**2)-2.9*alog(p)
c           bmx=sqrt(sig/10./pi)
c          endif
c         else
cc kmp
c          if(ish.ge.7)write(ifch,*)'sig_kmp chosen'
c          if(ecm.lt.2.5)then
c           call utindx(nkmp,kmpecm,ecm,iecm)
c           bmx=kmpbmx(iecm)
c          else
c           p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2)
c           sig=32.1+.66*(alog(p)**2)-5.6*alog(p)
c           bmx=sqrt(sig/10./pi)
c          endif
c         endif
c        elseif(kc(3).eq.-1)then    
cc strange particles involved
c         if(kc(1).eq.3.or.kc(2).eq.3)then
cc kpp
c          if(ish.ge.7)write(ifch,*)'sig_kpp chosen'
c          if(ecm.lt.2.5)then
c           call utindx(nkpp,kppecm,ecm,iecm)
c           bmx=kppbmx(iecm)
c          else
c           p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2)
c           sig=18.1+.26*(alog(p)**2)-alog(p)
c           bmx=sqrt(sig/10./pi)
c          endif
c         else
cc kpn
c          if(ish.ge.7)write(ifch,*)'sig_kpn chosen'
c           if(ecm.lt.2.5)then
c           call utindx(nkpn,kpnecm,ecm,iecm)
c           bmx=kpnbmx(iecm)
c          else
c           p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2)
c           sig=18.7+.21*(alog(p)**2)-.89*alog(p)
c           bmx=sqrt(sig/10./pi)
c          endif
c         endif
c        elseif(kc(3).ge.2)then         
cc two strange particles involved
c         bmx=1.
c        elseif(kc(3).le.-2)then 
cc two strange particles involved
c        bmx=1.
c        endif
c        
c       else !-----------baryon_cluster-anything----> 
c       
c        write(ifch,*)'i:',i,' id_i:',idptl(i),' j:',j,' id_j:',idptl(j)
c        write(ifch,*)'r_i:',radptl(i),' r_j:',radptl(j)
c        write(ifch,*)'E_cm:',ecm,' jc:',(jc(k,1),k=1,6),(jc(k,2),k=1,6)
c        write(ifch,*)'nq:',nq,' kc:',(kc(k),k=1,6)
c        bmx=radptl(i)+radptl(j)
c        
c       endif ! <--------------
c
cc -------------------------------------        
c      else  !   meson-meson
cc -------------------------------------        
c      
c       call idquad(i,nqi,nai,jci)
c       call idquad(j,nqj,naj,jcj)
c       do l=1,nflav
c        kci(l)=jci(l,1)-jci(l,2)
c        kcj(l)=jcj(l,1)-jcj(l,2)
c       enddo
c       
c       if(kci(3).eq.0.and.pptl(5,i).le.1.0
c     * .and.kcj(3).eq.0.and.pptl(5,j).le.1.0)then
c     
c        if(kci(1).eq.0.and.kci(2).eq.0.and.pptl(5,i).le.0.5)then
c         if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then
cc pi0 pi0
c          goto104
c         elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then
cc pi0 pi+
c          goto102
c         elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then
cc pi0 pi-
c          goto103
c         else
cc pi0 eta
c          goto105
c         endif
c        elseif(kci(1).eq.1.and.kci(2).eq.-1)then
c         if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then
cc pi+ pi0
c          goto102
c         elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then
cc pi+ pi+
c          goto101
c         elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then
cc pi+ pi-
c          goto104
c         else
cc pi+ eta
c          goto105
c         endif
c        elseif(kci(1).eq.-1.and.kci(2).eq.1)then
c         if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then
cc pi- pi0
c          goto103
c         elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then
cc pi- pi+
c          goto104
c         elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then
cc pi- pi-
c          goto101
c         else
cc pi- eta
c          goto105
c         endif
c        else
c         if(pptl(5,j).le.0.5)then
cc eta pi
c          goto105
c         else
cc eta eta
c          goto106
c         endif
c        endif
c101     continue
c        if(ish.ge.7)write(ifch,*)'sig_pi+pi+  chosen'
c        call utindx(npi1,pi1ecm,ecm,iecm)
c        bmx=pi1bmx(iecm)
c        goto110
c102     continue
c        if(ish.ge.7)write(ifch,*)'sig_pi+pi0  chosen'
c        call utindx(npi2,pi2ecm,ecm,iecm)
c        bmx=pi2bmx(iecm)
c        goto110
c103     continue
c        if(ish.ge.7)write(ifch,*)'sig_pi-pi0  chosen'
c        call utindx(npi3,pi3ecm,ecm,iecm)
c        bmx=pi3bmx(iecm)
c        goto110
c104     continue
c        if(ish.ge.7)write(ifch,*)'sig_pi-pi+  chosen'
c        call utindx(npi4,pi4ecm,ecm,iecm)
c        bmx=pi4bmx(iecm)
c        goto110
c105     continue
c        if(ish.ge.7)write(ifch,*)'sig_pi_eta  chosen'
c        call utindx(npi5,pi5ecm,ecm,iecm)
c        bmx=pi5bmx(iecm)
c        goto110
c106     continue
c        if(ish.ge.7)write(ifch,*)'sig_eta_eta  chosen'
c        bmx=0.4   ! approx.  sig=5mb
c110     continue
c
c       else !something else involved, strange etc.
c       
c        bmx=0.6  ! approx.  sig=10mb 
c        
c       endif
c       
cc --------------------------------         
c      endif
cc --------------------------------        
c
c      if(bij.le.bmx)ics=1
c      if(ish.ge.7)then
c      write(ifch,*)'b_mx:',bmx,' b_ij:',bij,' ics:',ics
c      endif
c
c1000  continue
c      call utprix('jintcs',ish,ishini,7)
c      return
c      end
c
cc----------------------------------------------------------------------
c      subroutine jintel(i,j,p,amim,xaver)
cc----------------------------------------------------------------------
cc  elastic scattering of ptls i,j
cc----------------------------------------------------------------------
c      include 'epos.inc'
c      real xaver(4)
c      real p(5),u(5),pei(5),pej(5)
c
cc     determine momenta of outgoing particles (pei,pej)
cc     -------------------------------------------------
c           if(p(5).le.(pptl(5,i)+pptl(5,j))*.99)then
c      if(ish.ge.1)then
c      call utmsg('jintel')
c      write(ifch,132)p(5),pptl(5,i)+pptl(5,j)
c132   format(1x,'*****  m_fus < m_i+m_j ---> qcm set zero    ( '
c     *,2f10.3,' )')
c      write(ifch,133)'p_i:  ',(pptl(k,i),k=1,5)
c      write(ifch,133)'p_j:  ',(pptl(k,j),k=1,5)
c      write(ifch,133)'p_fus:',(p(k),k=1,5)
c133   format(1x,a6,3x,5f10.3) 
c      call utmsgf
c      endif
c      qcm=0.
c           elseif(p(5).le.pptl(5,i)+pptl(5,j))then
c      qcm=0.
c           else
c      qcm=utpcm(p(5),pptl(5,i),pptl(5,j))
c           endif
c
cc isotropic
c      u(3)=2.*rangen()-1.
c      phi=2.*pi*rangen()
c      u(1)=sqrt(1.-u(3)**2)*cos(phi)
c      u(2)=sqrt(1.-u(3)**2)*sin(phi)
c      do 47 k=1,3
c      pei(k)= qcm*u(k)
c47    pej(k)=-qcm*u(k)
c
cc nonisotropic
cc-c   pt=sqrt(2./pi)*ptq*sqrt(-2*alog(rangen()) !2-dim Gauss
cc-c   if(pt.ge.qcm)pt=rangen()*qcm
cc-c   qpl=sqrt(qcm**2-pt**2)
cc-c   u(3)=qpl
cc-c   phi=2.*pi*rangen()
cc-c   u(1)=pt*cos(phi)
cc-c   u(2)=pt*sin(phi)
cc-c   call utaxis(i,j,a1,a2,a3)
cc-c   ivt=1
cc-c   if(a3.lt.0.)then
cc-c   a1=-a1
cc-c   a2=-a2
cc-c   a3=-a3
cc-c   ivt=-1
cc-c   endif 
cc-c   call utrota(-1,a1,a2,a3,u(1),u(2),u(3))
cc-c   do 47 k=1,3
cc-c   pei(k)= u(k)*ivt
cc-c47 pej(k)=-u(k)*ivt
c
c      pei(4)=sqrt(qcm**2+pptl(5,i)**2)
c      pej(4)=sqrt(qcm**2+pptl(5,j)**2)
c      pei(5)=pptl(5,i)
c      pej(5)=pptl(5,j)
c      call utlobo(-1,p(1),p(2),p(3),p(4),p(5)
c     *,pei(1),pei(2),pei(3),pei(4))
c      call utlobo(-1,p(1),p(2),p(3),p(4),p(5)
c     *,pej(1),pej(2),pej(3),pej(4))
c
cc     fill /cptl/
cc     -----------
c      do 49 lo=1,2
c      nptl=nptl+1
c      if(lo.eq.1)ij=i
c      if(lo.eq.2)ij=j
c      do 48 k=1,5
c      if(lo.eq.1)pptl(k,nptl)=pei(k)
c      if(lo.eq.2)pptl(k,nptl)=pej(k)
c48    continue
c      istptl(nptl)=0
c      idptl(nptl)=idptl(ij)
c      ibptl(1,nptl)=ibptl(1,ij)
c      ibptl(2,nptl)=ibptl(2,ij)
c      ibptl(3,nptl)=ibptl(3,ij)
c      ibptl(4,nptl)=ibptl(4,ij)
c      xorptl(1,nptl)=xaver(1)           
c      xorptl(2,nptl)=xaver(2)
c      xorptl(3,nptl)=xaver(3)
c      xorptl(4,nptl)=xaver(4)
c      iorptl(nptl)=i
c      jorptl(nptl)=j
c      tivptl(1,nptl)=xaver(4)                
c      call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
c      tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
c      ifrptl(1,nptl)=0
c      ifrptl(2,nptl)=0
c      ityptl(nptl)=ityptl(ij)
c49    continue
c
c1000  return
c      end
c
cc----------------------------------------------------------------------
c      subroutine jintfs(i,j,p,nq,jc,amim,iret) 
cc----------------------------------------------------------------------
cc  input:
cc    i,j: particle indices   
cc  output: 
cc    p: 5-momentum of fused ptl
cc    nq: net quark number of fused ptl
cc    jc: jc of fused ptl 
cc    amim : minimum mass of fused ptl
cc    iret: return code
cc            0=ok
cc            1=mass p(5) less than amim
cc----------------------------------------------------------------------
c      include 'epos.inc'
c      integer jci(nflav,2),jcj(nflav,2),jc(nflav,2)
c      real p(5)
c      double precision ppfus(4),pp52
c
c      iret=0
c
c      do 35 k=1,4
c       p(k)=pptl(k,i)+pptl(k,j)
c35    ppfus(k)=dble(p(k))
c      pp52=ppfus(4)**2-ppfus(3)**2-ppfus(2)**2-ppfus(1)**2
c      if(pp52.le.0.)then
c       if(ish.ge.1)then
c        call utmsg('jintfs')
c        write(ifch,*)'*****  mfus**2 < 0    (',pp52,' )'
c        write(ifch,*)(ppfus(m),m=1,4)
c        call utmsgf
c       endif
c       goto1001
c      endif
c      p(5)=sngl(dsqrt(pp52))
c
c      call idquad(i,nqi,nai,jci)
c      call idquad(j,nqj,naj,jcj)
c      do 29 n=1,nflav
c       do 29 k=1,2
c29    jc(n,k)=jci(n,k)+jcj(n,k)
c      nq=0
c      do 54 n=1,nflav
c54    nq=nq+jc(n,1)-jc(n,2)
c                       
c      call idcomj(jc)
c      amim=utamnz(jc,5)+amimfs
c      if(p(5).lt.amim)goto1001
c      goto1000
c1001  iret=1
c1000  return
c      end
c
cc----------------------------------------------------------------------
c      subroutine jintfu(i,j,p,jc)
cc----------------------------------------------------------------------
cc input:
cc   i,j: particle indices
cc   p,jc: momentum and jc of fused object
cc outout:
cc   id, ib, ist, ior, jor, ifr, ity of fused particle
cc   written to /cptl/ after increasing nptl by 1
cc----------------------------------------------------------------------
c
c      include 'epos.inc'
c      parameter (mxlook=10000,mxdky=2000)
c      common/dkytab/look(mxlook),cbr(mxdky),mode(5,mxdky)
c      real p(5)
c      integer jc(nflav,2),ic(2),ib(4)
c
c      nptl=nptl+1
c
cc momentum
c      do k=1,5
c       pptl(k,nptl)=p(k)
c      enddo
c      amf=p(5)
c
cc determine idr, ib(1-4)
c      idr=0
c      do 40 nf=1,nflav
c      do 40 ij=1,2
c      if(jc(nf,ij).ge.10)idr=7*10**8
c40    continue
c           if(idr/10**8.ne.7)then
c      call idenco(jc,ic,ireten)
c      if(ireten.eq.1)call utstop('jintfu: idenco ret code = 1&')
c      id=idtra(ic,0,0,3)
c43    amc=amf
c      call idres(id,amc,idr,iadj,1)
c          if(idr.ne.0)then
c      lid=look(iabs(idr))
c      if((lid.le.0.or.lid.gt.0.and.mode(2,lid).eq.0)
c     *.and.pptl(5,nptl).gt.amc+1e-3)then
c      amf=amf+0.010
c      goto43
c      endif
c      if((lid.le.0.or.lid.gt.0.and.mode(2,lid).eq.0)
c     *.and.abs(amc-pptl(5,nptl)).gt.1e-3)then
c      if(ish.ge.1)then
c      call utmsg('jintfu')
c      write(ifch,*)'*****  not on mass shell after fusion: '
c     *,pptl(5,nptl),amc
c      call utmsgf
c      endif
c      endif
c           endif
c           if(idr.eq.0)then
c      if(mod(ic(1),100).ne.0.or.mod(ic(2),100).ne.0)then
c      idr=9*10**8
c      else
c      idr=8*10**8+ic(1)*100+ic(2)/100
c      endif
c           endif
c           else
c      call idtrbi(jc,ib(1),ib(2),ib(3),ib(4))
c      idr=idr
c     *+mod(jc(1,1)+jc(2,1)+jc(3,1)+jc(4,1),10**4)*10**4
c     *+mod(jc(1,2)+jc(2,2)+jc(3,2)+jc(4,2),10**4)
c      if(ish.ge.7)write(ifch,*)'ib:',(ib(kk),kk=1,4)
c           endif
c
c      n=nptl
c
cc     fill /cptl/
cc     -----------
c      idptl(n)=idr
c      do k=1,4
c      ibptl(k,n)=ib(k)
c      enddo
c      istptl(n)=0
c      iorptl(n)=i
c      jorptl(n)=j
c      ifrptl(1,n)=0
c      ifrptl(2,n)=0
c      ityptl(n)=50
c
c      return
c      end
c
c----------------------------------------------------------------------
      subroutine jintpo
c----------------------------------------------------------------------
c  parton-ladder-fusion -- 3D grid -- based on string segments
c----------------------------------------------------------------------
      include 'epos.inc'
      include 'epos.incico'
      common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
     *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
      parameter(xgrid=8,mxgrid=10,kgrid=6,kegrid=3)
      parameter(ssgrid=0.5*xgrid*kgrid*kegrid)
      parameter(mssgrid=mxgrid*kgrid*kegrid)
      parameter(mxcl=1000,mxcli=200)
      integer idropgrid(mxgrid,mxgrid,mssgrid)
     &       ,jdropgrid(mxgrid,mxgrid,mssgrid)
     &       ,jclu(mssgrid),nsegp4(mxcl)
     &       ,irep(mxcl),mmji(mxcl,mxcli)
     &       ,jccl(mxcl,nflav,2),nseg(mxcl),mseg(mxcl),kclu(mxcl)
     &       ,naseg(0:mxcl),nfseg(mxcl),nsegmx(mxcl)
     &       ,ic(2),jc(nflav,2),ke(6),jcjj(nflav,2)
     &       ,nsegmt(mxptl)   !,nclk(mssgrid)
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
     &                ,pptld(5,mxcl),p4tmp
      common   /cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat 
      common/cdelzet/delzet,delsgr /cvocell/vocell
      common/cranphi/ranphi,ranecc,weiecc
      logical first
      double precision ptest(4),ttest,p52,p4mean(mxcl)  !,ppp(5),qqq(5),www(5)
      save ncntpo
      data ncntpo /0/
      ncntpo=ncntpo+1
      call utpri('jintpo',ish,ishini,4)

      if(ttaus.gt.kgrid)then
        write(ifmt,*)'ttaus,kgrid :',ttaus,kgrid
        stop'jintpo: ttaus too big.'
      endif
      tauico=ttaus
      fegrid=yhaha/5.36                 !rapidity range/rap range at RHIC
      sgrid=ttaus*ssgrid/kgrid *fegrid
      msgrid=max(1.4d0,ttaus)*mssgrid/kgrid *fegrid
      if(mod(msgrid,2).eq.0)msgrid=msgrid+1
      vocell= 2*xgrid/mxgrid * 2*xgrid/mxgrid * 2*sgrid/msgrid 
      delzet=2*sgrid/msgrid/ttaus
      delsgr=2*sgrid/msgrid
      nptla=nptl
      nsegsuj=max(3,nint(float(nsegsu)*fegrid))
      !print*,'dx,dz:',2*xgrid/mxgrid,2*sgrid/msgrid
      !. ,'  nx nz:',  mxgrid,msgrid, ttaus
                
c...compute x,y,z

      if(ish.ge.6)write(ifch,*)'compute x,y,z'
      do n=1,nptla
        if(istptl(n).eq.0)then
          call jtain(n,xptl(n),yptl(n),zptl(n),tptl(n),nnn,0)
          call jtaus(zptl(n),tz,sptl(n))
          strap=1e10
          xpl=tptl(n)+zptl(n)
          xmi=tptl(n)-zptl(n) 
          if(xmi.gt.0.0.and.xpl.gt.0.0)strap=0.5*log(xpl/xmi)    
          dezptl(n)=strap       !space-time-rapidity
c           if(iorptl(n).eq.4142)write(ifch,*)'ini',n
c     & ,pptl(4,n),pptl(3,n),idptl(n),xptl(n),yptl(n),zptl(n),sptl(n)
        else
          xptl(n)=0.
          yptl(n)=0.
          zptl(n)=0.
          sptl(n)=0.
          dezptl(n)=1e10
        endif
      enddo
      ntry=0

c...valid particles
 
      do n=1,nptla
        iaaptl(n)=1
        pt2=pptl(1,n)**2+pptl(2,n)**2
        am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2
        if(istptl(n).eq.0)then
          call idquac(n,idum1,idum2,idum3,jc)
        else
          jc(4,1)=0
          jc(4,2)=0
        endif
        if(istptl(n).ne.0.or.ityptl(n).ge.60)then
          iaaptl(n)=0        
	!elseif(ityptl(n).eq.41.or.ityptl(n).eq.51)then
	!  iaaptl(n)=0        !avoid particles coming from remnant droplet decay
	!   		      !(already droplet decay products !)
        elseif(abs(am2tmp-pptl(5,n)*pptl(5,n)).gt.30.)then
          iaaptl(n)=0           !to discard off shell particles
        elseif(jc(4,1).ne.0.or.jc(4,2).ne.0)then
          iaaptl(n)=0           !to discard charmed particles
        !elseif(pt2.lt.1.e-3)then
        !  iaaptl(n)=0           !to discard slow proton (spectators)
        elseif(pt2.gt.(1.5*ptclu)**2)then
          if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)iaaptl(n)=0
        elseif(pt2.gt.(0.5*ptclu)**2)then
          if(rangen().lt.(sqrt(pt2)-0.5*ptclu)/ptclu)then
            if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)iaaptl(n)=0
          endif
        endif
      enddo

c...Start cluster formation

 8888 continue
      ntry=ntry+1
      if(ntry.gt.90)        !do not put more than 100 or change limit for p4mean
     &call utstop('jintpo: cluster formation impossible ! &')
      nptl=nptla

      do 1 k=1,msgrid
      do 1 j=1,mxgrid
      do 1 i=1,mxgrid
      idropgrid(i,j,k)=0
  1   continue

c...count string segments in cell

      if(ish.ge.6)write(ifch,*)'count string segments in cell'
      do n=1,nptla
        if(iaaptl(n).ne.0)then
          i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid)
          j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid)
          k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid)
          if(  i.ge.1.and.i.le.mxgrid
     &    .and.j.ge.1.and.j.le.mxgrid
     &    .and.k.ge.1.and.k.le.msgrid)
     &      idropgrid(i,j,k)=idropgrid(i,j,k)+1
        endif
      enddo
  
c...check high pt segments
c      !to use this part one has to define:
c      !...  1 = valid
c      !... -1 = valid but high pt
c      !...  0 = not valid
c
c      if(ish.ge.6)write(ifch,*)'check high pt segments'
c      do n=1,nptla
c        if(iaaptl(n).eq.-1)then
c          i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid)
c          j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid)
c          k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid)
c          if(  i.ge.1.and.i.le.mxgrid
c     &    .and.j.ge.1.and.j.le.mxgrid
c     &    .and.k.ge.1.and.k.le.msgrid)then
c            if(idropgrid(i,j,k).lt.2*nsegce)then  !surface segments
c              iaaptl(n)=0
c              idropgrid(i,j,k)=idropgrid(i,j,k)-1
c              if(idropgrid(i,j,k).lt.0)stop'jintpo: not possible.    '
c           else
c             iaaptl(n)=1
c           endif
c          endif
c        endif
c      enddo
  
c...identify clusters

      if(ish.ge.6)write(ifch,*)'identify clusters'
      do k=1,msgrid !~~~~~~k-loop
        jjj=0
        do j=1,mxgrid
          first=.true.
          do i=1,mxgrid
            if(idropgrid(i,j,k).ge.nsegce)then
              if(first)then
                ifirst=i
                jjj=jjj+1
                if(jjj.gt.mxcl)stop'jintpo: mxcl too small.   ' 
                irep(jjj)=0
                jj=jjj
                first=.false.
              endif
              jdropgrid(i,j,k)=jj
              if(j.gt.1.and.jdropgrid(i,j-1,k).ne.0)then
                jjo=jdropgrid(i,j-1,k)
                if(jjo.lt.jj)then
                  if(jj.eq.jjj)jjj=jjj-1
                  jjx=jj
                  jj=jjo
                  do ii=ifirst,i
                    jdropgrid(ii,j,k)=jj
                    if(jdropgrid(ii,j-1,k).eq.jjx)then
                      if(jjx.gt.jjj)jjj=jjj+1
                        jja=jjx
                      jjb=jj
  90                  continue 
                      if(irep(jja).eq.0.or.irep(jja).eq.jjb)then
                        irep(jja)=jjb
                      else
                        mn=min(irep(jja),jjb)
                        mx=max(irep(jja),jjb)
                        irep(jja)=mn
                        jja=mx
                        jjb=mn
                        goto90
                      endif
                    endif
                  enddo
                elseif(jdropgrid(i,j-1,k).gt.jj)then
                  irep(jjo)=jj
                endif  
              endif
            else
              jdropgrid(i,j,k)=0
              first=.true.
            endif  
          enddo
        enddo
        !~~~~cluster jj ---> cluster irep(jj)       
        do jj=jjj,1,-1
         if(irep(jj).ne.0)then
           do j=1,mxgrid
             do i=1,mxgrid
               if(jdropgrid(i,j,k).eq.jj)jdropgrid(i,j,k)=irep(jj)
             enddo
           enddo
         endif
        enddo
        !~~~~~remove empty cluster indices       
        jjjx=jjj   
        jjj=0
        jj=0
        do jjx=1,jjjx
         if(irep(jjx).eq.0)then
           jj=jj+1
           jjj=jjj+1
         else
           do j=1,mxgrid
             do i=1,mxgrid
               if(jdropgrid(i,j,k).gt.jj)
     &           jdropgrid(i,j,k)=jdropgrid(i,j,k)-1
             enddo
           enddo
         endif
        enddo        
        !~~~~~        
        jclu(k)=jjj      
      enddo !~~~~~~~~~~~~~~~~~ END k-loop

c...absolute clusters numbering (for all k)

      if(ish.ge.6)write(ifch,*)'absolute clusters numbering'
      jjj=jclu(1)
      do k=2,msgrid
        do j=1,mxgrid
          do i=1,mxgrid
            if(jdropgrid(i,j,k).gt.0)
     &        jdropgrid(i,j,k)=jdropgrid(i,j,k)+jjj
          enddo
        enddo  
        jjj=jjj+jclu(k)
      enddo
      jjs=0
      do k=1,msgrid
        do jj=1,jclu(k)
          jjs=jjs+1
          kclu(jjs)=k
         nseg(jjs)=0
         nsegp4(jjs)=0
          p4mean(jjs)=0d0
        enddo
      enddo
      nptlb=nptl+jjj

c...calculate mean energy of segments going into in clusters
         
      if(ish.ge.6)write(ifch,*)
     &'calculate mean energy of segments going into in clusters'
      do 95 n=1,nptla
        if(iaaptl(n).eq.0)goto 95
        i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid)
        j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid)
        k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid)
        if(    i.ge.1.and.i.le.mxgrid
     &          .and.j.ge.1.and.j.le.mxgrid
     &          .and.k.ge.1.and.k.le.msgrid)then
          jj=jdropgrid(i,j,k)
          if(jj.gt.0)then
            nsegp4(jj)=nsegp4(jj)+1
            p4mean(jj)=p4mean(jj)+dble(pptl(4,n))
c           if(jj.eq.9)write(ifch,*)'after',n
c     & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k
         endif  
        endif   
 95   continue

      do jj=1,jjj
        if(nsegp4(jj).gt.5)then
          p4mean(jj)=p4mean(jj)/dble(nsegp4(jj))
        else
          p4mean(jj)=1d50
        endif
      enddo
      
c...mark segments going into in clusters, count them
         
      if(ish.ge.6)write(ifch,*)'mark segments going into in clusters'
      do 96 n=1,nptla
        if(iaaptl(n).eq.0)goto96
        i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid)
        j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid)
        k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid)
        if(    i.ge.1.and.i.le.mxgrid
     &          .and.j.ge.1.and.j.le.mxgrid
     &          .and.k.ge.1.and.k.le.msgrid)then
          jj=jdropgrid(i,j,k)
          if(jj.gt.0)then
          ! count only particles with reasonnable energy
            pt2=pptl(1,n)**2+pptl(2,n)**2+pptl(5,n)**2
            am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2
            if(abs(am2tmp).lt.0.1.or.
     &         dble(pptl(4,n)).le.100.d0/dble(ntry)*p4mean(jj))then
              istptl(n) = 3 
              nseg(jj)=nseg(jj)+1
c           if(jj.eq.9)write(ifch,*)'after',n
c     & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k
            else
c              write(ifmt,*)'Rejected particle : ',n,idptl(n)
             if(ish.ge.1)write(ifch,*)'Rejected particle : ',n,idptl(n)
     & ,ityptl(n),(pptl(k,n),k=1,5),am2tmp
     & ,dble(pptl(4,n))/dble(nsegp4(jj)),p4mean(jj),nsegp4(jj),jj
              nsegp4(jj)=nsegp4(jj)-1
              iaaptl(n)=0        
              if(nsegp4(jj).gt.5)then
                p4mean(jj)=(p4mean(jj)*dble(nsegp4(jj)+1)-pptl(4,n))
     &               /dble(nsegp4(jj))
              else
                p4mean(jj)=1d50
              endif
            endif
         endif  
        endif   
  96  continue


c...add segments moving towards clusters
         
      if(ish.ge.6)write(ifch,*)'add segments moving towards clusters'
      if(iocluin.eq.1)then
      do 93 n=1,nptla
        if(iaaptl(n).eq.0)goto93
        ihit=0
        i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid)
        j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid)
        k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid)
        if(    i.ge.1.and.i.le.mxgrid
     &          .and.j.ge.1.and.j.le.mxgrid
     &          .and.k.ge.1.and.k.le.msgrid)then
          jgr=jdropgrid(i,j,k)
          if(jgr.eq.0)then
           if(i.ge.mxgrid/2)then
            isi=-1
           else
            isi=1
           endif   
           if(j.ge.mxgrid/2)then
            jsi=-1
           else
            jsi=1
           endif   
            do ii=i,i+2*isi,isi
             do jj=j,j+2*jsi,jsi
               if(.not.(ii.eq.i.and.jj.eq.j))then
                if(ii.ge.1.and.ii.le.mxgrid
     .          .and.jj.ge.1.and.jj.le.mxgrid)then
                 jg=jdropgrid(ii,jj,k)
                 if(jg.gt.0)then
                  if(nseg(jg).gt.50)then
                   x=xptl(n)
                   y=yptl(n)
                   vrad=( x*pptl(1,n)/pptl(4,n)+y*pptl(2,n)/pptl(4,n)) 
                    if(vrad.lt.0.)then
                     ihit=1
                     goto94
                    endif
                  endif
                 endif
                endif
               endif 
             enddo
            enddo
          endif  
        endif
   94   continue        
        if(ihit.eq.1)then
          istptl(n) = 3 
         nseg(jg)=nseg(jg)+1
         delx=2*xgrid/mxgrid*(ii-i)
         dely=2*xgrid/mxgrid*(jj-j)
         xptl(n)=xptl(n)+delx
         yptl(n)=yptl(n)+dely
        ! if(k.ge.10.and.k.le.13)
   !  .          print*,k, '  ',i,j,'  ',ii,jj,nseg(jg),'   ',delx,dely
        endif   
  93  continue
      endif
      
c...prepare /cptl/ for clusters

      if(ish.ge.6)write(ifch,*)'prepare /cptl/ for clusters'
      do jj=1,jjj
         nptl=nptl+1
          istptl(nptl)=12
          do l=1,4
            pptl(l,nptl)=0.
            xorptl(l,nptl)=0
          enddo
          sptl(nptl)=0
          uptl(nptl)=0
          optl(nptl)=0
          desptl(nptl)=0
c limit the maximum number of subcluster to half the number of particle
c (not to have empty subclusters)
          nsegmx(jj)=max(1,nint(float(nseg(jj))/float(nsegsuj)))
      enddo

c...prepare /cptl/ for subclusters

      if(ish.ge.6)write(ifch,*)'prepare /cptl/ for subclusters'
      mm=0
      do jj=1,jjj
        do ii=1,nsegmx(jj)
          mm=mm+1
          if(mm.gt.mxcl)stop'jintpo: mxcl too small.        '
          mmji(jj,ii)=mm
          mseg(mm)=0
          nptl=nptl+1
          istptl(nptl)=10
          do l=1,4
            pptld(l,mm)=0d0
            pptl(l,nptl)=0.
            xorptl(l,nptl)=0
          enddo
          sptl(nptl)=0
          uptl(nptl)=0
          optl(nptl)=0
          desptl(nptl)=0
          do l=1,nflav
            jccl(mm,l,1)=0
            jccl(mm,l,2)=0
          enddo
          iorptl(nptl)=nptla+jj
          jorptl(n)=0 
          if(ii.eq.1)ifrptl(1,nptla+jj)=nptlb+mm
          ifrptl(2,nptla+jj)=nptlb+mm
        enddo  
      enddo

c...separate string segments, add dense area segments to clusters

      if(ish.ge.6)write(ifch,*)'separate string segments'
      do 98 n=1,nptla
        if(istptl(n).ne.3)goto 98
        i=1+(xptl(n)+xgrid)/(2*xgrid/mxgrid)
        j=1+(yptl(n)+xgrid)/(2*xgrid/mxgrid)
        k=1+(sptl(n)+sgrid)/(2*sgrid/msgrid)
        if(    i.ge.1.and.i.le.mxgrid
     &          .and.j.ge.1.and.j.le.mxgrid
     &          .and.k.ge.1.and.k.le.msgrid)then
          jj=jdropgrid(i,j,k)
          if(jj.gt.0)then
            iimx=nsegmx(jj)
            do ii=1,iimx
              mm=mmji(jj,ii)
              if(mseg(mm).eq.0)goto 10          !not to have an empty cluster
            enddo
            ii=1+rangen()*iimx
            ii=min(ii,iimx)
            iini=ii
            am2tmpmx=1e30
 9          ntmp=mmji(jj,ii)
            am2tmp=(pptld(4,ntmp)+pptld(3,ntmp))
     &            *(pptld(4,ntmp)-pptld(3,ntmp))
            if(am2tmp.gt.5.e3)then
              if(am2tmp.lt.am2tmpmx)then
                mm=mmji(jj,ii)
                am2tmpmx=am2tmp
              endif
              ii=ii+1
              if(ii.gt.iimx)ii=1
              if(ii.ne.iini)then
                goto 9
              else
                goto 10
              endif
            endif
            mm=mmji(jj,ii)
 10         continue
            mseg(mm)=mseg(mm)+1
            ifrptl(1,n)=mm      !local use of ifrptl
c           write(ifch,*)'end',n
c     & ,pptl(4,n),pptl(3,n),idptl(n),i,j,k,sptl(n)
            p4tmp=0d0
            do l=1,3
             pptld(l,mm)=  pptld(l,mm)  + dble(pptl(l,n))
             p4tmp=p4tmp+dble(pptl(l,n))*dble(pptl(l,n))
            enddo
            p4tmp=sqrt(p4tmp+dble(pptl(5,n)*pptl(5,n)))
            pptld(4,mm)=  pptld(4,mm)  + p4tmp
c           if(mm.eq.86)write(ifch,*)'other',n
c     & ,pptl(4,n),pptl(3,n),pptl(5,n),idptl(n),k,p4tmp
c     & ,pptld(4,mm),pptld(3,mm),pptld(2,mm),pptld(1,mm)
c     & ,(pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm))
            xorptl(1,nptlb+mm)=xorptl(1,nptlb+mm)+xptl(n)
            xorptl(2,nptlb+mm)=xorptl(2,nptlb+mm)+yptl(n)
            xorptl(3,nptlb+mm)=xorptl(3,nptlb+mm)+zptl(n)
            xorptl(4,nptlb+mm)=xorptl(4,nptlb+mm)+tptl(n)
            sptl(nptlb+mm)=sptl(nptlb+mm)+sptl(n)
            aa=cos(phievt)
            bb=sin(phievt)
            cc=-sin(phievt)
            dd=cos(phievt)
            xrot=xptl(n)*aa+yptl(n)*bb
            yrot=xptl(n)*cc+yptl(n)*dd
            uptl(nptlb+mm)=uptl(nptlb+mm)+xrot**2
            optl(nptlb+mm)=optl(nptlb+mm)+yrot**2
            desptl(nptlb+mm)=desptl(nptlb+mm)+xrot*yrot
            id=idptl(n)
            ida=iabs(id/10)
            ids=id/iabs(id)
            if(ida.ne.111.and.ida.ne.222.and.ida.ne.333)id=id/10*10
            if(ida.eq.111.or. ida.eq.222.or. ida.eq.333)id=id/10*10+ids
            if(ida.eq.213)id=1230*ids
            ic(1)=idtrai(1,id,1)
            ic(2)=idtrai(2,id,1)
            call iddeco(ic,jc)
            do l=1,nflav
              jccl(mm,l,1)=jccl(mm,l,1)+jc(l,1)
              jccl(mm,l,2)=jccl(mm,l,2)+jc(l,2)
            enddo
          else
            idropgrid(i,j,k)=0
          endif  
        endif
  98  continue
  
c...associate segments to clusters

      if(ish.ge.6)write(ifch,*)'associate segments to clusters'
      naseg(0)=0
      do jj=1,jjj
        do ii=1,nsegmx(jj)
          mm=mmji(jj,ii)
          naseg(mm)=naseg(mm-1)+mseg(mm)
          nfseg(mm)=0
        enddo  
      enddo
      do 97 n=1,nptla
        if(istptl(n).ne.3)goto97
        istptl(n) = 2
        mm=ifrptl(1,n)
        nfseg(mm)=nfseg(mm)+1
        nsegmt(naseg(mm-1)+nfseg(mm))=n
 97   enddo        
      do jj=1,jjj
        nst=0
        do ii=1,nsegmx(jj)
          mm=mmji(jj,ii)
          if(mseg(mm).ne.nfseg(mm))stop'jintpo: mseg.ne.nfseg        '
          nst=nst+mseg(mm)
        enddo
        if(nst.ne.nseg(jj))stop'sum(mseg(mm)).ne.nseg(jj)'
      enddo

c...finish cluster storage to /cptl/

      if(ish.ge.6)write(ifch,*)'finish cluster storage to /cptl/'
      xx=0.
      yy=0.
      xy=0.
      mjjsegsum=0
      do jj=1,jjj
       njj=nptla+jj 
       mjjseg=0
       do l=1,nflav
         jcjj(l,1)=0
         jcjj(l,2)=0
       enddo
       do ii=1,nsegmx(jj)
        mm=mmji(jj,ii)
        n=nptlb+mm

        do l=1,nflav
          jc(l,1)=jccl(mm,l,1)
          jc(l,2)=jccl(mm,l,2)
          ke(l)=jc(l,1)-jc(l,2)
          jcjj(l,1)=jcjj(l,1)+jc(l,1)
          jcjj(l,2)=jcjj(l,2)+jc(l,2)
        enddo
        call idenct(jc,idptl(n)
     *  ,ibptl(1,n),ibptl(2,n),ibptl(3,n),ibptl(4,n))
        ttest=0d0
       do ji=1,4
         ptest(ji)=0d0
          do ns=naseg(mm-1)+1,naseg(mm)
            ni=nsegmt(ns)
            ptest(ji)=ptest(ji)+dble(pptl(ji,ni))
          enddo
         ptest(ji)=abs(ptest(ji)-pptld(ji,mm))
         ttest=ttest+ptest(ji)
        enddo
        p52=(pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm))
     &      -pptld(1,mm)**2-pptld(2,mm)**2 
        pptld(5,mm)=0d0
        if(p52.gt.0d0)then
         pptld(5,mm)=sqrt(p52)
         jerr(2)=jerr(2)+1
       elseif(p52.le.0d0)then
        jerr(3)=jerr(3)+1
        pptld(5,mm)=dble(utamnu(ke(1),ke(2),ke(3),ke(4),ke(5),ke(6),5))
         pptld(4,mm)=sqrt(pptld(3,mm)*pptld(3,mm)
     &                   +pptld(2,mm)*pptld(2,mm)
     &                   +pptld(1,mm)*pptld(1,mm)
     &                   +pptld(5,mm)*pptld(5,mm))
        endif
      if(ish.ge.1.and.(abs(ttest).gt.1.d0.or.pptld(5,mm).gt.120d0))then
         call utmsg('jintpo&')
         write(ifmt,*)'***** Warning in jintpo !',ntry
         write(ifch,*)'***** jintpo: momenta messed up (ttest > 0)'
          write(ifch,*)'*****',mm,n,mseg(mm),p52,ttest
          write(ifch,*)'*****',jj,nsegp4(jj),p4mean(jj)
         write(ifch,'(a,10x,4f15.4)')'*****',(pptld(ji,mm),ji=1,4)
          do ns=naseg(mm-1)+1,naseg(mm)
            ni=nsegmt(ns)
            write(ifch,'(a,2i5,5f15.4,f12.4)')'*****',ni,idptl(ni)
     *     ,(pptl(ji,ni),ji=1,4),pptl(5,ni)**2
     *     ,(pptl(4,ni)+pptl(3,ni))*(pptl(4,ni)-pptl(3,ni))
     *       -pptl(1,ni)**2-pptl(2,ni)**2 
          enddo
       endif
       if(pptld(5,mm).gt.120d0)then
          p4max=0.
          nh=0
          do ns=naseg(mm-1)+1,naseg(mm)
            ni=nsegmt(ns)
            if(pptl(4,ni).ge.p4max)then
              nh=ni
              p4max=pptl(4,ni)
            endif
          enddo
          if(nh.le.0)then
            stop'Cannot be in jintpo ...'
          else   !put back nh as normal particle 
            iaaptl(nh)=0        
            istptl(nh)=0        
          endif
         if(ish.ge.1)
     &    write(ifch,*)'***** Redo cluster without heavy particle :'
     &                 ,nh,ntry
          goto 8888
       endif
        do l=1,5
         pptl(l,n)=sngl(pptld(l,mm))
        enddo
       mjjseg=mjjseg+mseg(mm)
        do l=1,4
          pptl(l,njj)=pptl(l,njj)+pptl(l,n)
          xorptl(l,njj)=xorptl(l,njj)+xorptl(l,n)
          xorptl(l,n)=xorptl(l,n)/float(mseg(mm))
        enddo
        sptl(njj)=sptl(njj)+sptl(n)
        uptl(njj)=uptl(njj)+uptl(n)
        optl(njj)=optl(njj)+optl(n)
        desptl(njj)=desptl(njj)+desptl(n)
        sptl(n)=sptl(n)/float(mseg(mm))
        uptl(n)=uptl(n)/float(mseg(mm))
        optl(n)=optl(n)/float(mseg(mm))
        desptl(n)=desptl(n)/float(mseg(mm))
        radptl(n)=0
        istptl(n)=10
        ifrptl(1,n)=0 
        ifrptl(2,n)=0 
        tivptl(1,n)=xorptl(4,n)
        tivptl(2,n)=ainfin
        ityptl(n)=60
        vocri=pptl(5,n)/epscri(ioclude)
        vo=max(vocri,vocell)
        radptl(n)=(3.*vo/4./pi)**0.3333
        dezptl(n)=0.
       enddo
       do l=1,4
        xorptl(l,njj)=xorptl(l,njj)/float(mjjseg)
       enddo
       mjjsegsum=mjjsegsum+mjjseg
       xx=xx+uptl(njj)
       yy=yy+optl(njj)
       xy=xy+desptl(njj)
       sptl(njj)=sptl(njj)/float(mjjseg)
       uptl(njj)=uptl(njj)/float(mjjseg)
       optl(njj)=optl(njj)/float(mjjseg)
       desptl(njj)=desptl(njj)/float(mjjseg)
       pjj52=(pptl(4,njj)+pptl(3,njj))*(pptl(4,njj)-pptl(3,njj))
     &    -pptl(1,njj)**2-pptl(2,njj)**2
        pptl(5,njj)=0
        if(pjj52.gt.0)pptl(5,njj)=sqrt(pjj52)
        ityptl(njj)=60
        call idenct(jc,idptl(njj)
     *  ,ibptl(1,njj),ibptl(2,njj),ibptl(3,njj),ibptl(4,njj))
         !.................
         !inertia tensor:   
         !<y**2>   -<x*y>  with <x**2>=uptl  <y**2>=optl
         !-<x*y>   <x**2>                   <xy>=desptl
         !.................
      enddo
      if(mjjsegsum.gt.0)then
        xx=xx/float(mjjsegsum)
        yy=yy/float(mjjsegsum)
        xy=xy/float(mjjsegsum)
        dta=0.5*abs(xx-yy)
        ev1=(xx+yy)/2+sqrt(dta**2+xy**2)
        ev2=(xx+yy)/2-sqrt(dta**2+xy**2)
        if(xy.lt.0..and.dta.ne.0.)then
          ranphi=0.5*atan(-xy/dta)
        elseif(xy.gt.0..and.dta.ne.0.)then
          ranphi=-0.5*atan(xy/dta)
        else
          ranphi=0
        endif
      else
        ranphi=0
      endif
                 
c...print
      
      if(ish.ge.5)then
        write(ifch,*)'print'
        do k=1,msgrid
        write(ifch,*)'k=',k,'  jclu=',jclu(k)
        do j=mxgrid,1,-1
        write(ifch,'(10i4,3x,10i4)')(idropgrid(i,j,k),i=1,mxgrid)
     &                ,(jdropgrid(i,j,k),i=1,mxgrid)    
        enddo
        enddo
        write(ifch,'(a,a)')
     &    '    k   jj  nseg      mm  mseg     n      mass'
     &   ,'         s         y         z         t '
        do jj=1,jjj
         do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
           mm=mmji(jj,ii)
           n=nptlb+mm
           sg=pptl(3,n)/abs(pptl(3,n))
           tm=sqrt(pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2)
           y=sg*alog((pptl(4,n)+sg*pptl(3,n))/tm)
c           if(kclu(jj).eq.44)print *,tm,pptl(4,n),pptl(3,n),iorptl(n)
           write(ifch,'(2i5,i6,i8,2i6,5f10.3)')
     &       kclu(jj),jj,nseg(jj),mm,mseg(mm),n,pptl(5,n)
     &      ,sptl(n),y,xorptl(3,n),xorptl(4,n)
         enddo
        enddo
      endif
             
c...compute ranecc

      weiecc=0
      ranecc=0
      do jj=1,jjj
       do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
        mm=mmji(jj,ii)
        ip=nptlb+mm
          ipo=iorptl(ip)
          xx=uptl(ipo)   ! <x**2>
          yy=optl(ipo)   ! <y**2>
          xy=desptl(ipo) ! <x*y>
          dta=0.5*abs(xx-yy)
          ev1=(xx+yy)/2+sqrt(dta**2+xy**2)
          ev2=(xx+yy)/2-sqrt(dta**2+xy**2)
          yy=ev1
          xx=ev2
              ecc=(yy-xx)/(yy+xx)
          ranecc=ranecc+ecc*pptl(5,ip)
          weiecc=weiecc+pptl(5,ip)
       enddo
      enddo
      if(weiecc.gt.0.)then
       ranecc=ranecc/weiecc
      endif

c...decay
             
      if(ish.ge.6)write(ifch,*)'decay'
      if(ifrade.eq.0.or.ispherio.gt.0)goto1000
      if(jdecay.eq.0)goto1000
ctp060829      npl=0
      do jj=1,jjj
       do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
        mm=mmji(jj,ii)
        np=nptlb+mm
        nptlij=nptl
        iret=1
c        ishb=ish
c        ish=0
        call hnbaaa(np,iret)
c        iret=1
c       ish=ishb
        if(iret.eq.1)then
          istptl(np)=istptl(np)+2
          do ns=naseg(mm-1)+1,naseg(mm)
            n=nsegmt(ns)
            istptl(n) = 0  
          enddo
        else 
          istptl(np)=istptl(np)+1
          ifrptl(1,np)=nptlij+1
          ifrptl(2,np)=nptl
          do n=nptlij+1,nptl
            iorptl(n)=np
            jorptl(n)=0
            istptl(n)=0
            ifrptl(1,n)=0
            ifrptl(2,n)=0
              !xorptl(1,n)  already defined in hnbaaa
              !xorptl(2,n)  already defined in hnbaaa
              !xorptl(3,n)  already defined in hnbaaa
              !xorptl(4,n)  already defined in hnbaaa
            tivptl(1,n)=xorptl(4,n)
            call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
            r=rangen()
            tivptl(2,n)=tivptl(1,n)+taugm*(-alog(r))
              !ityptl(n)   already defined in hnbaaa
            radptl(n)=0.
            dezptl(n)=0.
            itsptl(n)=0
            rinptl(n)=kclu(jj)-msgrid/2
          enddo
        endif  
       enddo
      enddo
       
 1000 continue           
      call utprix('jintpo',ish,ishini,4)
      end

cc-----------------------------------------------------------------------
c      subroutine jrad(i,nq,na,jc,rad)
cc-----------------------------------------------------------------------
cc     return hadron radius (data taken from huefner and povh)
cc-----------------------------------------------------------------------
c      include 'epos.inc'
c      integer jc(nflav,2),kc(nflav)
c
c      id=iabs(idptl(i))
c      am=pptl(5,i)
c      if(id.lt.10000)then
c       k=mod(id,10)
c      else
c       k=1
c      endif
c      do l=1,nflav
c       kc(l)=iabs(jc(l,1)-jc(l,2))
c      enddo
c
c      if(nq.eq.0)then   ! mesons
c       if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0.and.kc(4).eq.0)then
c        if(k.eq.0)then           ! flavor singlet pseudoscalar mesons
c         if(am.ge.0.000)then
c          rad=0.64                 ! pi0
c          if(am.ge.0.500)then
c           rad=0.60                ! eta
c           if(am.ge.0.900)then
c            rad=0.40               ! eta prime
c            if(am.ge.2.900)then
c             rad=0.17              ! eta charm
c            endif
c           endif
c          endif
c         else
c          write(ifch,*)
c     *    'i:',i,' id:',idptl(i),' k:',k,' m:',am
c          write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
c          call utstop('jrad: meson radius not defined&')
c         endif 
c        else                    ! flavor singlet vector mesons
c         if(am.ge.0.000)then
c          rad=0.72                ! rho,omega
c          if(am.ge.1.000)then
c           rad=0.46               ! phi
c           if(am.ge.3.000)then
c            rad=0.20              ! J/psi
c           endif
c          endif
c         else
c          write(ifch,*)
c     *    'i:',i,' id:',idptl(i),' k:',k,' m:',am
c          write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
c          call utstop('jrad: meson radius not defined&')
c         endif
c        endif
c       elseif(kc(3).eq.0.and.kc(4).eq.0)then  ! nonstrange, noncharmed
c        if(k.eq.0)then
c         rad=0.64  ! pi
c        else
c         rad=0.72  ! resonances
c        endif
c       elseif(kc(3).ne.0.and.kc(4).eq.0)then  ! strange
c        if(k.eq.0)then
c         rad=0.59  ! kaons
c        else
c         rad=0.68  ! kaon resonances
c        endif
c       else                                   ! charmed
c        write(ifch,*)'i:',i,' id:',idptl(i)
c        call utstop('jrad: radius of meson not defined&')
c       endif
c      else   !baryons
c       if(kc(4).gt.0)then       ! charmed
c        write(ifch,*)
c     *  'i:',i,' id:',idptl(i),' k:',k,' m:',am
c        write(ifch,*)'i:',i,' id:',idptl(i)
c        call utstop('jrad: radius of charmed baryon not defined&')
c       elseif(kc(3).eq.0)then   ! nonstrange
c        if(k.eq.0)then
c         rad=0.82  !nucleons
c        else
c         rad=1.00  !resonances
c        endif
c       elseif(kc(3).eq.1)then   ! strange
c        if(k.eq.0)then
c         rad=0.76  !lambda, sigma
c        else
c         rad=0.93  !resonances
c        endif
c       elseif(kc(3).eq.2)then   ! double strange
c        if(k.eq.0)then
c         rad=0.71  !cascades
c        else
c         rad=0.87  !resonances
c        endif
c       elseif(kc(3).ge.3)then   ! triple strange
c        rad=0.79  !omega
c       else
c        write(ifch,*)
c     *  'i:',i,' id:',idptl(i),' k:',k,' m:',am
c        write(ifch,*)
c     *  'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6),
c     *  ' |q-qbar|:',(kc(l),l=1,6)
c        call utstop('jrad: should not happen&')
c       endif
cc string fragments with |#q|>3
c       if(na.gt.3)then
c        a=(na/3.)**(1./3.)
c        if(ish.ge.7)then
c         call utmsg('jrad ')
c         write(ifch,*)
c     *   'i:',i,' id:',idptl(i),' k:',k,' m:',am
c         write(ifch,*)
c     *   'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6),
c     *   ' |q-qbar|:',(kc(l),l=1,6)
c         write(ifch,*)'nq:',nq,' na:',na,' r:',rad,' ar:',a*rad
c         call utmsgf
c        endif
c        rad=rad*a
c       endif
c      endif
c
c      if(ish.ge.7)then
c       write(ifch,*)
c     * 'i:',i,' id:',idptl(i),' k:',k,' m:',am,' rad:',rad
c       write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
c      endif
c
c      return
c      end
c
c-----------------------------------------------------------------------
      subroutine jresc
c-----------------------------------------------------------------------
      include 'epos.inc'
      double precision pa(5),pj(5)
      integer ipptl(mxptl)

      call utpri('jresc ',ish,ishini,4)

      iret=0
      nptlpt=maproj+matarg
      np=0
      do i=nptlpt+1,nptl
       if(istptl(i).eq.0
     * .and.idptl(i).lt.10000.and.pptl(5,i).gt.0.01)then
        np=np+1
        ipptl(np)=i
       endif
      enddo
      if(np.lt.2)goto1001
      do ii=1,np
       i=ipptl(ii)
       if(mod(iabs(idptl(i)),10).lt.2)then
        call idmass(idptl(i),ami)
        dm=abs(ami-pptl(5,i))
        if(dm.gt.0.001)then
         ntry=0
1        continue
         ntry=ntry+1
2        jj=1+int(rangen()*np)
         j=ipptl(jj)
         if(ish.ge.4)write(ifch,*)i,j,istptl(j)
         if(j.eq.i)goto2
         if(mod(iabs(idptl(j)),10).lt.2)then
          call idmass(idptl(j),amj)
         else
          amj=pptl(5,j)
         endif
         do l=1,5
          pa(l)=dble(pptl(l,i))
          pj(l)=dble(pptl(l,j))
         enddo
         if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70)
         if(ish.ge.4)write(ifch,11)i,idptl(i),'before:',pa,'want:',ami
         if(ish.ge.4)write(ifch,11)j,idptl(j),'before:',pj,'want:',amj
         call jrescl(pa,dble(ami),pj,dble(amj),iret)
         if(iret.eq.1)then
          if(ntry.le.50)then
           goto1
          else
           goto1001
          endif
         endif
         if(ish.ge.4)write(ifch,11)i,idptl(i),' after:',pa
         if(ish.ge.4)write(ifch,11)j,idptl(j),' after:',pj
         if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70)
         do l=1,5
          pptl(l,i)=sngl(pa(l))
          pptl(l,j)=sngl(pj(l))
         enddo
        endif
       endif
      enddo
11    format(i5,1x,i5,1x,a,1x,5(d8.2,1x),a,1x,e8.2)

1000  continue 
      call utprix('jresc ',ish,ishini,4)
      return

1001  continue 
      if(ish.ge.1)then
        write(ifmt,'(a)')'jresc: could not put on shell'
      endif
      goto1000

      end

c-----------------------------------------------------------------------
      subroutine jrescl(p1,am1,p2,am2,iret)
c-----------------------------------------------------------------------
c rescale momenta of two particles such that the masses assume given 
c values.
c input:
c   p1, p2: momenta of the two particles
c   am1, am2: desired masses of the two particles
c output:
c   p1, p2: new momenta of the two particles
c-----------------------------------------------------------------------
      include 'epos.inc'
      double precision p1(5),p2(5)
     *                ,p1n(5),p2n(5)
     *                ,a1,a2,a12,am1,am2
     *                ,b1,b2,c,d,e,f,g,p,q,r

      call utpri('jrescl',ish,ishini,7)

      iret=0
      a1=p1(5)**2
      a2=p2(5)**2
      a12=p1(4)*p2(4)-p1(3)*p2(3)-p1(2)*p2(2)-p1(1)*p2(1)
      if(a12.le.(a1+a2))then
       if(ish.ge.7)write(ifch,*)'a_12 < a_1 + a_2'
       if(ish.ge.7)write(ifch,*)a12,' < ',a1+a2
c      goto1001
      endif

11    format(5(d9.3,1x))
      if(ish.ge.7)write(ifch,11)p1,a1
      if(ish.ge.7)write(ifch,11)p2,a2
      if(ish.ge.7)write(ifch,*)a12

      c=(a1+a12)/(a2+a12)
      d=(a1-am1**2-a2+am2**2)/(a2+a12)*0.5d0

      e=a1-2d0*a12*c+a2*c**2
      f=2d0*(a1-a12*(c+d)+a2*c*d)
      g=a1-2d0*a12*d+a2*d**2-am1**2

      p=f/e
      q=g/e
      r=p**2-4d0*q

      if(ish.ge.7)write(ifch,*)'c:',c,' d:',d
      if(ish.ge.7)write(ifch,*)'e:',e,' f:',f,' g:',g
      if(ish.ge.7)write(ifch,*)'p:',p,' q:',q,' r:',r
      if(r.lt.0d0)goto1001

      b1=-0.5d0*(p-dsqrt(r))

      b2=b1*c+d

      if(ish.ge.7)write(ifch,*)'b_1:',b1,' b_2:',b2

      do i=1,4
       p1n(i)=(1d0+b1)*p1(i)-b2*p2(i)
       p2n(i)=(1d0+b2)*p2(i)-b1*p1(i)
      enddo

      a1=p1n(4)**2-p1n(3)**2-p1n(2)**2-p1n(1)**2
      a2=p2n(4)**2-p2n(3)**2-p2n(2)**2-p2n(1)**2
      if(a1.gt.0d0.and.a2.gt.0d0)then
       do i=1,4
        p1(i)=p1n(i)
        p2(i)=p2n(i)
       enddo
       p1(5)=dsqrt(a1)
       p2(5)=dsqrt(a2)
       if(ish.ge.7)write(ifch,11)p1,a1
       if(ish.ge.7)write(ifch,11)p2,a2
      else
       goto1001
      endif

      if(p1(4).lt.0..or.p2(4).lt.0.)goto1001

1000  continue 
      call utprix('jrescl',ish,ishini,7)
      return

1001  continue 
      iret=1
      goto1000
      end

c-----------------------------------------------------------------------
      subroutine jtain(i,x,y,z,t,n,iopt)
c-----------------------------------------------------------------------
c returns intersection (x,y,z,t) of ptl-i-trajectory with taus-line.
c input:    
c   i: particle number
c   iopt: formation time considered (0) or not (1) 
c output:
c   x,y,z,t: 4-vector of intersection point
c   n: exit code
c       n=0: ok
c       n=1: ptl lives later
c       n=2: ptl lives earlier
c       n=9: tiv1>tiv2
c-----------------------------------------------------------------------
      include 'epos.inc'
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat 
      double precision vv,zza,zz,tt,xo3,xo4,ti1,ti2,derr,dd
      double precision ttp,zzp,ttt,zzt,vvt,vvp,spt2m2E,p4
      common/ctfi/tin,tfi
      double precision ttau0
      common/cttau0/ttau0
      
      n=0

      tin=0
      tfi=0

      derr=1d-2
      ttp=tpro*ttaus
      zzp=zpro*ttaus
      ttt=ttar*ttaus
      zzt=ztar*ttaus
      vv=sign(min(1.d0,abs(dble(pptl(3,i)))/dble(pptl(4,i)))
     &                ,dble(pptl(3,i)))


      if(abs(vv).ge.1.d0)then
        spt2m2E=dble(pptl(1,i)*pptl(1,i)+pptl(2,i)*pptl(2,i)
     &              +pptl(5,i)*pptl(5,i))
c        if(pptl(4,i).le.0.)then
          p4=sqrt(dble(pptl(3,i)*pptl(3,i))+spt2m2E)
c        else
c          p4=dble(pptl(4,i))
c        endif
ctp to avoid precision problem, replace abs(p3)/p4 by sqrt(1-(pt2+m2)/E2)
        spt2m2E=min(1.d0,sqrt(spt2m2E)/p4)
        vv=sign(sqrt((1d0+spt2m2E)*(1d0-spt2m2E)),dble(pptl(3,i)))
      endif
      xo3=dble(xorptl(3,i))
      xo4=dble(xorptl(4,i))
      zza=xo3-xo4*vv
      if(iopt.eq.0)ti1=dble(tivptl(1,i))
      if(iopt.eq.1)ti1=dble(xo4)
      ti2=dble(tivptl(2,i))
      
      if(ti1.gt.ti2)then
        n=9
        goto1
      endif
      
      zfi=sngl(xo3+(ti2-xo4)*vv)
      call jtaus(zfi,tzfi,szfi)
      tfi=tzfi
      if(tfi.ge.sngl(ti2))then
        n=2
        goto1
      endif
 
      zin=sngl(xo3+(ti1-xo4)*vv)
      call jtaus(zin,tzin,szin)
      tin=tzin
      if(tin.le.sngl(ti1))then
        n=1
        goto1
      endif


    1 continue
 
           if(ttaus.le.ttau0)then
      tt=ttaus
      zz=xo3+(tt-xo4)*vv
      if(tt.lt.ti1.and.n.eq.0)n=1
      if(tt.ge.ti2.and.n.eq.0)n=2
      goto1000
           else
      vvt=zzt/ttt
      vvp=zzp/ttp
      tt=(ttt+(zza-zzt)*vvt)/(1-vv*vvt)
      zz=xo3+(tt-xo4)*vv
      if(zz.le.zzt)then
      if(tt.lt.ti1.and.n.eq.0)n=1
      if(tt.ge.ti2.and.n.eq.0)n=2
      goto1000
      endif
      tt=(ttp+(zza-zzp)*vvp)/(1-vv*vvp)
      zz=xo3+(tt-xo4)*vv
      if(zz.ge.zzp)then
      if(tt.lt.ti1.and.n.eq.0)n=1
      if(tt.ge.ti2.and.n.eq.0)n=2
      goto1000
      endif
      dd=1-vv**2
      if(sngl(dd).eq.0..and.vv.gt.0.)then
      tt=-(ttaus**2+zza**2)/2d0/zza
      elseif(sngl(dd).eq.0..and.vv.lt.0.)then
      tt=(ttaus**2+zza**2)/2d0/zza
      else
      tt=(zza*vv+dsqrt(zza**2+ttaus**2*dd))/dd
      endif
      zz=xo3+(tt-xo4)*vv
      if(tt.lt.ti1.and.n.eq.0)n=1
      if(tt.ge.ti2.and.n.eq.0)n=2
        if(dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr*ttaus**2.and.
     *dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr)then
      if(ish.ge.1)then
      call utmsg('jtain')
      write(ifch,*)'*****  ttaus**2 .ne. (tt+zz)*(tt-zz)'
      write(ifch,*)sngl(ttaus**2),sngl((tt+zz)*(tt-zz))
      call utmsgf
      endif
      goto1000
        endif
           endif
 
1000  t=sngl(tt)
      z=sngl(zz)
      x=xorptl(1,i)+(t-xorptl(4,i))*pptl(1,i)/pptl(4,i)
      y=xorptl(2,i)+(t-xorptl(4,i))*pptl(2,i)/pptl(4,i)
      return
      end

c-----------------------------------------------------------------------
      subroutine jtaix(i,tau,zor,tor,z,t)
c-----------------------------------------------------------------------
c     returns intersection z,t of ptl-i-trajectory with hyperbola h.
c        h: (t-tor)**2-(z-zor)**2=tau**2 .
c        zor, tor double precision.
c-----------------------------------------------------------------------
      include 'epos.inc'
      double precision tor,zor,tors,zors,vv,cc,dd,ttau,derr,tt,zz
      derr=1d-3
      ttau=dble(tau)
      zors=dble(xorptl(3,i))-zor
      tors=dble(xorptl(4,i))-tor
      vv=dble(pptl(3,i)/pptl(4,i))
      vv=dmin1(vv,1d0)
      vv=dmax1(vv,-1d0)
      cc=zors-tors*vv
      dd=1d0-vv**2
      dd=dmax1(dd,0d0)
           if(dd.eq.0d0.and.cc.eq.0d0)then
      if(tau.eq.0.)tt=0d0
      if(tau.ne.0.)tt=dble(ainfin)
      zz=tt
      goto1000
           elseif(dd.eq.0d0)then
      tt=-(ttau**2+cc**2)/2d0/cc/vv
           elseif(dd.lt.1e-8)then
      tt=-(ttau**2+cc**2)/2d0/cc/vv
      call utmsg('jtaix')
      write(ifch,*)'*****  dd = ',dd,'    treated as zero'
      call utmsgf
           else
      tt=(cc*vv+dsqrt(cc**2+ttau**2*dd))
      tt=tt/dd
           endif
      zz=cc+tt*vv
      if(dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr*ttau**2
     *.and.dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr
     *.and.tors**2-zors**2.lt.1e6)then
      if(ish.ge.1)then
      call utmsg('jtaix')
      write(ifch,*)'*****  ttau**2 .ne. (tt+zz)*(tt-zz)'
      write(ifch,*)sngl(ttau**2),sngl((tt+zz)*(tt-zz))
      write(ifch,*)'tau,t,z:'
      write(ifch,*)tau,tt,zz
      write(ifch,*)'#,id(ptl):',i,idptl(i)
      write(ifch,*)'zor,tor(str):',zor,tor
      write(ifch,*)'zors,tors,p,e(ptl):'
      write(ifch,*)sngl(zors),sngl(tors),pptl(3,i),pptl(4,i)
      call utmsgf
      endif
      endif
1000  z=sngl(zz+zor)
      t=sngl(tt+tor)
      return
      end

c-----------------------------------------------------------------------
      subroutine jtaug(su,so,g,y)
c-----------------------------------------------------------------------
c  returns g factor and rapidity y for given su, so
c-----------------------------------------------------------------------
      include 'epos.inc'
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat 
      double precision ttp,zzp,ttt,zzt,ssp,sst,ssu,sso,ss1,ss2,gg
     *,ssav,yyav,hh
      double precision ttau0
      common/cttau0/ttau0

      ssu=dble(su)
      sso=dble(so)

      if(ssu.ge.sso)then 
      sso=(ssu+sso)*0.5d0 + dble(abs(dezzer))*ttaus*0.5d0
      ssu=(ssu+sso)*0.5d0 - dble(abs(dezzer))*ttaus*0.5d0
      so=real(sso)
      su=real(ssu)
      endif
      if(ssu.ge.sso)then
        print*,ssu,sso,dble(abs(dezzer))*ttaus*0.5d0
        stop'STOP: sr jtaug: ssu.ge.sso'
      endif  
      
      g=1

      if(ttaus.le.ttau0)return

      ttp=tpro*ttaus    
      zzp=zpro*ttaus
      ttt=ttar*ttaus
      zzt=ztar*ttaus
      ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp))
      sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt))

      ssav=(ssu+sso)/2d0
      yyav=ssav/ttaus
      if(ssav.ge.ssp)yyav=detap 
      if(ssav.le.sst)yyav=detat

      gg=0
      if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu)
      if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu))
      if(ssu.lt.ssp.and.sso.gt.sst)then
      ss1=dmax1(ssu,sst)
      ss2=dmin1(sso,ssp)
      gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) )
      endif
      gg=gg/(sso-ssu)

      hh=0
      if(ssu.lt.sst)hh=hh + dsinh(detat-yyav) * (dmin1(sst,sso)-ssu)
      if(sso.gt.ssp)hh=hh + dsinh(detap-yyav) * (sso-dmax1(ssp,ssu))
      if(ssu.lt.ssp.and.sso.gt.sst)then
      ss1=dmax1(ssu,sst)
      ss2=dmin1(sso,ssp)
      hh=hh+ttaus*( dcosh(ss2/ttaus-yyav)-dcosh(ss1/ttaus-yyav) )
      endif
      hh=hh/(sso-ssu)

      yyav=yyav+0.5d0*dlog((gg+hh)/(gg-hh))

      gg=0
      if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu)
      if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu))
      if(ssu.lt.ssp.and.sso.gt.sst)then
      ss1=dmax1(ssu,sst)
      ss2=dmin1(sso,ssp)
      gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) )
      endif
      gg=gg/(sso-ssu)

      g=sngl(gg)
      y=sngl(yyav)
      return
      end

c-----------------------------------------------------------------------
      subroutine jtaui(s,ts,zs)
c-----------------------------------------------------------------------
c  returns time ts and coord zs corresponding to ttaus and inv. length s
c-----------------------------------------------------------------------

      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat 
      double precision ttau0
      common/cttau0/ttau0

      double precision ttp,zzp,ttt,zzt,ssp,sst,ss,zeta

      zs=s
      ts=sngl(ttaus)

      if(ttaus.le.ttau0)return

      ttp=tpro*ttaus    
      zzp=zpro*ttaus
      ttt=ttar*ttaus
      zzt=ztar*ttaus
      ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp))
      sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt))
      ss=dble(s)

           if(ss.le.sst)then
      zs=sngl(zzt+ttar*(ss-sst))
      ts=sngl(ttt+(dble(zs)-zzt)*zzt/ttt)
           elseif(ss.ge.ssp)then
      zs=sngl(zzp+tpro*(ss-ssp))
      ts=sngl(ttp+(dble(zs)-zzp)*zzp/ttp)
           else
      zeta=ss/ttaus
      ts=sngl(ttaus*dcosh(zeta))
      zs=sngl(ttaus*dsinh(zeta))
           endif
           
      return
      end

c-----------------------------------------------------------------------
      subroutine jtauin
c-----------------------------------------------------------------------
c initializes equal time hyperbola at ttaus
c-----------------------------------------------------------------------
      include 'epos.inc'
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat 
      double precision ttau0,rcproj,rctarg
      common/geom1/rcproj,rctarg
      common/cttau0/ttau0

      call utpri('jtauin',ish,ishini,6)

      if(ttaus.gt.ttau0)then
       if(rcproj.gt.1d-10)then
        detap=dmin1(dble((ypjtl-yhaha)*etafac),dlog(ttaus/rcproj))
       else
        detap=dble((ypjtl-yhaha)*etafac)
       endif
       if(rctarg.gt.1d-10)then
        detat=dmax1(dble(-yhaha*etafac),dlog(rctarg/ttaus))
       else
        detat=dble(-yhaha*etafac)
       endif
       tpro=dcosh(detap)
       zpro=dsinh(detap)
       ttar=dcosh(detat)
       ztar=dsinh(detat)
      else
       detap=0d0
       detat=0d0
       tpro=0d0
       zpro=0d0
       ttar=0d0
       ztar=0d0
      endif

      if(ish.ge.6)then
       write(ifch,*)'hyperbola at tau=',ttaus
       write(ifch,*)'r_p:',rcproj,' r_t:',rctarg
       write(ifch,*)'y_p:',detap,' y_t:',detat
       write(ifch,*)'t_p:',tpro,' z_p:',zpro
       write(ifch,*)'t_t:',ttar,' z_t:',ztar
      endif

      call utprix('jtauin',ish,ishini,6)
      return
      end

c-----------------------------------------------------------------------
      subroutine jtaus(z,tz,sz)
c-----------------------------------------------------------------------
c  returns time tz and inv length sz corresponding to ttaus and z
c-----------------------------------------------------------------------
      include 'epos.inc'
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat 
      double precision ttau0
      common/cttau0/ttau0

      double precision ttp,zzp,ttt,zzt,zz,tzz

      tz=sngl(ttaus)
      sz=z

      if(ttaus.le.ttau0)return

      ttp=tpro*ttaus
      zzp=zpro*ttaus
      ttt=ttar*ttaus
      zzt=ztar*ttaus
      zz=dble(z)

           if(zz.le.zzt)then
      tz=sngl(ttt+(zz-zzt)*zzt/ttt)
      sz=sngl(ttaus*detat+(zz-zzt)/ttar)
           elseif(zz.ge.zzp)then
      tz=sngl(ttp+(zz-zzp)*zzp/ttp)
      sz=sngl(ttaus*detap+(zz-zzp)/tpro)
           else
      if(sngl(ttaus).ge.ainfin)then
      tz=sngl(ttaus)
      sz=0.
      if(ish.ge.1)then
      call utmsg('jtaus')
      write(ifch,*)'*****  large ttaus; set tz=ttaus, sz=0'
      write(ifch,*)'ttaus=',ttaus,'zz=',zz
      call utmsgf
      endif
      else
      tzz=dsqrt(ttaus**2+zz**2)
      tz=sngl(tzz)
      sz=sngl(ttaus*0.5d0*dlog((tzz+zz)/(tzz-zz)))
      endif
           endif

      return
      end

c-----------------------------------------------------------------------
      subroutine jtaux(t,z,ttaux)
c-----------------------------------------------------------------------
c  returns ttaux (-> tau-line) for given t and z.
c  ttaux: double precision.
c-----------------------------------------------------------------------
      double precision ttaux,tt,zz,rcproj,rctarg,zt1,zp1,zt2,zp2,ttau0
      common/geom1/rcproj,rctarg
      common/cttau0/ttau0
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/   tpro,zpro,ttar,ztar,ttaus,detap,detat 

      tt=dble(t)
      zz=dble(z)

      if(tt.gt.ttau0)then
       zt1=rctarg-tt
       zp1=tt-rcproj
       zt2=ztar/ttar*tt
       zp2=zpro/tpro*tt
       if(zz.le.dmax1(zt1,zt2))then
        if(zt1.gt.zt2)then
         ttaux=rctarg*dsqrt((tt-zz)/(2d0*rctarg-tt-zz))
        else
         ttaux=(ttar*tt-ztar*zz)/(ttar**2-ztar**2)
        endif
       elseif(zz.ge.dmin1(zp1,zp2))then
        if(zp1.lt.zp2)then
         ttaux=rcproj*dsqrt((tt+zz)/(2d0*rcproj-tt+zz))
        else
         ttaux=(tpro*tt-zpro*zz)/(tpro**2-zpro**2)
        endif
       else
        ttaux=dsqrt(tt**2-zz**2)
       endif
      else
       ttaux=tt
      endif

      return
      end

c-----------------------------------------------------------------------
      subroutine xjden1(ii,itau,x,y,rad,o,u)
c-----------------------------------------------------------------------
c ii=0: initialization
c ii=1: determining dense regions in space of individual events
c       x,y,rad: tranverse coordinates and radius of particle i
c       o,u: specifies long range: u < s < o (s: long coordinate)
c ii=2: plot of individual event
c       xpar1: itau ; valid: 1,..,10
c       xpar2: z-range: -xpar2 < z < xpar2 
c       xpar3, x-range: -xpar3 < x < xpar3 
c       xpar4, y-range: -xpar4 < y < xpar4 
c-----------------------------------------------------------------------
      include "epos.inc"
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/   tpro,zpro,ttar,ztar,ttaus,detap,detat 

      if(idensi.ne.1)stop'STOP in xjden1: idensi must be set 1'
    
      dlcoox=0.5
      dlcooy=0.5

           if(ii.eq.0)then

      do i=1,nzeta
      do j=1,mxcoox
      do k=1,mxcooy
      kdensi(itau,i,j,k)=0
      enddo
      enddo
      enddo

           elseif(ii.eq.1)then

      if(itau.lt.1.or.itau.gt.mxtau)return

      tau=sngl(ttaus)
      zu=u/tau
      zo=o/tau

            do 1 i=1,nzeta
      zi=-nzeta*dlzeta/2-dlzeta/2+i*dlzeta
      if(zu.gt.zi.or.zo.lt.zi)goto1
      do 2 j=1,mxcoox
      xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox
      do 3 k=1,mxcooy
      yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
      if((x-xj)**2+(y-yk)**2.gt.rad**2)goto3
      kdensi(itau,i,j,k)=1
    3 continue
    2 continue
    1 continue

           elseif(ii.eq.2)then

      itaux=nint(xpar1)
      if(itaux.gt.mxtau)stop'STOP in xjden1: itaux too large'

      iz=nint(xpar2/dlzeta)
      ix=nint(xpar3/dlcoox)
      iy=nint(xpar4/dlcooy)
      if(iz.gt.nzeta/2)stop'STOP in xjden1: zeta-range too large'
      if(ix.gt.mxcoox/2)stop'STOP in xjden1: x-range too large'
      if(iy.gt.mxcooy/2)stop'STOP in xjden1: y-range too large'

      do k=mxcooy/2+1-iy,mxcooy/2+iy
      write(ifhi,'(a,f7.2)')      '! tau: ',tauv(itaux)
      write(ifhi,'(a)')      'openhisto'
      write(ifhi,'(a,2f7.2)')'xrange',-iz*dlzeta,iz*dlzeta
      write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox
      write(ifhi,'(a)')      'set ityp2d 3'
      write(ifhi,'(a)') 'txt  "xaxis space-time rapidity [z]"'
      write(ifhi,'(a)') 'txt  "yaxis transverse coordinate x (fm)"'
      write(ifhi,'(a,i4)')   'array2d',2*iz
      do j=mxcoox/2+1-ix,mxcoox/2+ix
      write(ifhi,'(40i2)')    (kdensi(itaux,i,j,k),
     *                        i=nzeta/2+1-iz,nzeta/2+iz) 
      enddo
      write(ifhi,'(a)')       '  endarray'
      write(ifhi,'(a)')       'closehisto plot2d'
      enddo

           else
      
      stop'STOP in xjden1: wrong option'

           endif

      return
      end

c-----------------------------------------------------------------------
      subroutine xjden2(ii,itau,x,y,rad,s)
c-----------------------------------------------------------------------
c ii=0: initialization
c ii=1: determining dense regions in space of individual events
c       x,y,rad: tranverse coordinates and radius of particle i
c       s: long coordinate
c ii=2: plot of individual event
c       xpar1: itau ; valid: 1,..,10
c       xpar2: s-range: -xpar2 < s < xpar2 
c       xpar3, x-range: -xpar3 < x < xpar3 
c       xpar4, y-range: -xpar4 < y < xpar4 
c-----------------------------------------------------------------------
      include "epos.inc"
      double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      common/cttaus/   tpro,zpro,ttar,ztar,ttaus,detap,detat 
      parameter (mxcoos=60)
      common/cdensh/kdensh(matau,mxcoos,mxcoox,mxcooy),ktot(matau)
      character cy*3

      dlcoox=0.5
      dlcooy=0.5
      dlcoos=0.5

           if(ii.eq.0)then

      do i=1,mxcoos
      do j=1,mxcoox
      do k=1,mxcooy
      kdensh(itau,i,j,k)=0
      enddo
      enddo
      enddo
      ktot(itau)=0

           elseif(ii.eq.1)then

      if(itau.lt.1.or.itau.gt.mxtau)return

      tau=sngl(ttaus)
      z=s/tau

      do 1 i=1,mxcoos
      si=-mxcoos*dlcoos/2-dlcoos/2+i*dlcoos
      do 2 j=1,mxcoox
      xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox
      do 3 k=1,mxcooy
      yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
      if(((x-xj)**2+(y-yk)**2+(z-si)**2).gt.rad**2)goto3
      kdensh(itau,i,j,k)=kdensh(itau,i,j,k)+1
      ktot(itau)=ktot(itau)+1
    3 continue
    2 continue
    1 continue

           elseif(ii.eq.2)then

      itaux=nint(xpar1)
      if(itaux.gt.mxtau)stop'STOP in xjden2: itaux too large'

      is=nint(xpar2/dlcoos)
      ix=nint(xpar3/dlcoox)
      iy=nint(xpar4/dlcooy)
      if(is.gt.mxcoos/2)stop'STOP in xjden2: s-range too large'
      if(ix.gt.mxcoox/2)stop'STOP in xjden2: x-range too large'
      if(iy.gt.mxcooy/2)stop'STOP in xjden2: y-range too large'

      do k=mxcooy/2+1-iy,mxcooy/2+iy
      write(cy,'(f3.1)')-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
      write(ifhi,'(a)')      'openhisto'
      write(ifhi,'(a,2f7.2)')'xrange',-is*dlcoos,is*dlcoos
      write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox
      write(ifhi,'(a)')      'set ityp2d 5'
      write(ifhi,'(a)') 'txt  "xaxis [z] "'
      write(ifhi,'(a)')
     *'txt  "yaxis  x (fm), y='//cy//' fm"'
      write(ifhi,'(a,i4)')   'array2d',2*is
      do j=mxcoox/2+1-ix,mxcoox/2+ix
      do i=mxcoos/2+1-is,mxcoos/2+is
      write(ifhi,'(e11.3)')
     *kdensh(itaux,i,j,k)/dlcooy/dlcoos/dlcoox/ktot(itaux)
      enddo
      enddo
      write(ifhi,'(a)')       '  endarray'
      write(ifhi,'(a)')       'closehisto plot2d'
      enddo

           else
      
      stop'STOP in xjden2: wrong option'

           endif

      return
      end

cc-----------------------------------------------------------------------
c      subroutine postscript(iii,ii,ic)
cc-----------------------------------------------------------------------
c      include 'epos.inc'
c      character*10 color(10)
c      if(iii.eq.0)then
c      ifps=21
c      open(unit=ifps,file='zzz.ps',status='unknown')
c      WRITE(ifps,'(a)') '%!PS-Adobe-2.0'
c      WRITE(ifps,'(a)') '%%Title: tt2.fig'
c      WRITE(ifps,'(a)') '%%Orientation: Portrait'
c      WRITE(ifps,'(a)') '%%BeginSetup'
c      WRITE(ifps,'(a)') '%%IncludeFeature: *PageSize A4'
c      WRITE(ifps,'(a)') '%%EndSetup'
c      WRITE(ifps,'(a)') '%%EndComments'
c      WRITE(ifps,*) '/l {lineto} bind def'
c      WRITE(ifps,*) '/rl {rlineto} bind def'
c      WRITE(ifps,*) '/m {moveto} bind def'
c      WRITE(ifps,*) '/rm {rmoveto} bind def'
c      WRITE(ifps,*) '/s {stroke} bind def'
c      WRITE(ifps,*) '/gr {grestore} bind def'
c      WRITE(ifps,*) '/gs {gsave} bind def'
c      WRITE(ifps,*) '/cp {closepath} bind def'
c      WRITE(ifps,*) '/tr {translate} bind def'
c      WRITE(ifps,*) '/sc {scale} bind def'
c      WRITE(ifps,*) '/sd {setdash} bind def'
c      WRITE(ifps,*) '/sdo {[.01 .05] 0 sd} bind def'
c      WRITE(ifps,*) '/sdf {[1 .0] 0 sd} bind def'
c      WRITE(ifps,*) '/n {newpath} bind def'
c      WRITE(ifps,*) '/slw {setlinewidth } bind def'
c      write(ifps,*) '/srgb {setrgbcolor} bind def'
c      write(ifps,*) '/lgrey      { 0 0.95 0.95 srgb} bind def'
c      write(ifps,*) '/black      { 0 0 0 srgb} bind def'
c      write(ifps,*) '/red        { 1 0 0 srgb} bind def  '
c      write(ifps,*) '/green      { 0 1 0  srgb} bind def  '
c      write(ifps,*) '/blue       { 0 0 1  srgb} bind def  '
c      write(ifps,*) '/yellow     { 1 0.5 0  srgb} bind def  '
c      write(ifps,*) '/turquoise  { 0 1 1  srgb} bind def  '
c      write(ifps,*) '/purple     { 1 0 1  srgb} bind def  '
cc      write(ifps,*) '/  {   srgb} bind def  '
cc      write(ifps,*) '/  {   srgb} bind def  '
c      write(ifps,*) '/ef {eofill} bind def'
c      WRITE(ifps,'(a)') '%%EndProlog'
c      WRITE(ifps,*) 'gsave'
c      WRITE(ifps,*) '/Helvetica findfont 10 scalefont setfont'
c      color(9)='lgrey     '
c      color(1)='black     '
c      color(2)='red       '
c      color(3)='green     '
c      color(4)='blue      '
c      color(7)='yellow    '
c      color(5)='turquoise '
c      color(6)='purple    '
c      np=0
c         elseif(iii.eq.1)then
c      np=np+1
c      write(ifps,'(a,i4)') '%%Page: number ',np
c      write(ifps,'(a)') 'gsave'
c      WRITE(ifps,*) '100 700 tr'
c      scale=0.125
c      WRITE(ifps,*) 1./scale,1./scale,' sc'
c      WRITE(ifps,*) scale/2.,' slw'
c      WRITE(ifps,*) '/Helvetica findfont ',15.*scale
c     &     ,' scalefont setfont'
c      write(ifps,*) color(1),' n ',smin,xmin,' m ( tau:',tau,') show '
c
c      WRITE(ifps,*) '/Helvetica findfont ',5.*scale
c     &     ,' scalefont setfont'
c
c
c      yb=-2.
c      dy=4./12.
c      yb=yb-dy/2
c      do iyb=0,11
c        yb=yb+dy
c        WRITE(ifps,*) 'gsave'
c        WRITE(ifps,*) (xmax-xmin)*1.1*float(int(iyb/4))
c     &       ,-(xmax-xmin)*1.1*mod(iyb,4),' tr'
c        write(ifps,*) ' n ',smin,xmin,' m ',smax,xmin,' l '
c     &       ,smax,xmax,' l ',smin,xmax,' l cp s '
cc.......particles in layer iyb.............
c        do i=1,nptl
c          if(ii.gt.0)then
c              write(ifps,*)  color(mod(i,5)+2)
c     &             ,' n ',u,x-r,' m ',o,x-r,' l '
c     &             ,o,x+r,' l ',u,x+r,' l cp s '
c              write(ifps,*) ' n ',u,x-r,' m (',i,ior,') show ' 
c          else
c              write(ifps,*) ' n ',s,x,r,0,360,' arc ',color(ic),' s '
c              write(ifps,*) ' n ',s-r,x,' m (',i,io,') show ' 
c          endif
c 10     enddo
c        write(ifps,*) color(1),' n ',smin,xmin,' m (',yb,') show ' 
c        WRITE(ifps,*) 'grestore'
c      enddo                    !yb bin
c      write(ifps,'(a)') 'grestore'
c      write(ifps,*) 'showpage'        
c          elseif(iii.eq.2)then
c       write(ifps,*) 'gr' 
c      
c       write(ifps,'(a)') '%%Trailer'
c       write(ifps,'(a,i4)') '%%Pages: ',np
c       write(ifps,'(a)') '%%EOF'
c       close(unit=ifps)
c          endif
c
c      return
c      end
c

c----------------------------------------------------------------------------------------
      subroutine xtauev(iii)
c----------------------------------------------------------------------------------------   
      end
c----------------------------------------------------------------------------------------
      subroutine wimi
c----------------------------------------------------------------------------------------   
      end
c----------------------------------------------------------------------------------------
      subroutine wimino
c----------------------------------------------------------------------------------------   
      end
c----------------------------------------------------------------------------------------
      subroutine xspace(iii)
c----------------------------------------------------------------------------------------   
      end
c----------------------------------------------------------------------------------------
      subroutine wclu
c----------------------------------------------------------------------------------------   
      end
c----------------------------------------------------------------------------------------
      subroutine wclufi
c----------------------------------------------------------------------------------------   
      end
c----------------------------------------------------------------------------------------
      subroutine wtime(iii)
c----------------------------------------------------------------------------------------   
      end