      Subroutine tpicm (kappa)
c *** see RH Landau, Program lpott, 1981
c     Subroutine tpicm calcs  on-shell pin  c.m. amplitude
c     kappa= pi-nucleon c.m. momentum
c     intialize If phase shifts  change
      Implicit Real*8 (a-h,i,k,m,o-z)
      Dimension nifty(20), ecm(100), d(14,100), eta(14,100), kap(100), t
     1(14), tx(14), bn(16), bnuc(20), rewave(14,100), imwave(14,100),
     2 bnf(16), bnucf(20)
      Common /sec2/ bnuc,bn,bnucf,bnf,t,tx,hbarc,pi,mpi,mn,nz,nes,nwaves
     1,nifty,na
      Common /nlspfl/ gin(14,100),kapg(100),denom(14,25),ecmden(25),fmn,
     1rewave,imwave,kap,nkapgs,necms
ccc   Equivalence (rewave(1,1),d(1,1)), (imwave(1,1),eta(1,1))
c     new fit to low energy phase shifts,using martin salomon s fit up t
c     and rhl fit to d and f waves (eta==1,ok till ca 150mev/c)
c     a*s are dummy variables,proprtnal to alpha (.5(exp(2idelta)-1))/2
c     code for bn, bnf, bnuc, bnucf
c     all for pi plus, reaarange for pi minus latter
c     1 res,p 2 imsp    3 res,n        4 imsn
c     5 rep,p 6imp,p    4 rep,n        8 imp,n
c     9 red,p 10imd,p   11 red,n       12 imd,n
c     13 ref,p          14 imf,p       15 ref,n  16 imf,n
c     17 reg,p          18 img,p       19 reg,n  20 img,n
c     b(pi+p)=sum(j) (j+1/2))t3/2(j)
c     b(pi+n)=sum(j) (j+1/2)*(t3/2(j)/3. + 2*r1/2(j)/3)
c     flip
c     bj(pi+p)=t3/2(j=l+1/2)- t3/2(-)
c     bf(pi+n)=1/3 (t3/2(+) -t3/2(-)) +2/3 (t1/2(+) -t1/2(-)
c     spin flip amplitudes bnf have the same code(1-4==0)
c     t(pi+p)= t(pi-n)  , t(pi+n)=t(pi-p)
c     twave(nwave,ne) = contrbtn to partl wve t matrx element from wave
c     at pi-n c.m. energy kap(ne)
c        nwave code
c     1-----s31     2-----p31     3-----p33     4-----s11
c     5-----p11     6-----p13     7-----d33     8-----d35
c     9-----f35     10----f37     11----d13     12----d15
c     13----f15     14----f17
c
      data ndata/1/
      If (kappa.le.-100.) ndata = 1
      kappa = abs(kappa)
      If (ndata.eq.0) GoTo 110
      kappa = abs(kappa)
      If (nifty(8).ne.2.and.nifty(8).ne.3) GoTo 100
c     nifty(8) .eq. 2,or 3 Read in phase shifts,cals t*s,store
c     If folding,scattering lengths/ms used for very low e
c     tmatrix(unfolded) for this nuceles,t folderd latter
c     amplitudes from input phase shifts and enrgy
c     ecm in mev, delta (d) Read in in degrees, f calculated in mev-1
c     eta=exp(-2 imdelta)  , If Read in as 0 then set to 1.
c
      nei = 1
ccc   If (nifty(8).ne.1.and.nifty(8).ne.3) GoTo 30
      If (nifty(8).ne.3) GoTo 30
c     nifty(8) =1or 3 1,set up low e delta   needed for folding
      nei = 26
      nes = nes+25
      Read (1,170) ecm(1)
      Rewind 1
      s = ecm(1)**2
      kappa = Sqrt((s-(mn+mpi)*(mn+mpi))*(s-(mn-mpi)*(mn-mpi))/4./s)
      If (kappa.le.5.) kappa = 6.
      kaph = (kappa-5.)/20.
      kappa = 0.
      Do 20 nkap=1,25
         kappa = kappa+kaph
         If (nkap.le.6) kappa = (nkap*1.5)
         Do 10 ni=1,14
            Call tpirsl (kappa,ni,ddum,w,tr,ti)
            eta(ni,nkap) = 1.
            d(ni,nkap) = ddum*180./pi
            t(ni) = tr
            tx(ni) = ti
 10      Continue
         ecm(nkap) = w
 20   Continue
 30   Write (6,180) nes,nwaves
      Read (1,170) (ecm(ne),(eta(nwave,ne),d(nwave,ne),nwave=1,nwaves),
     1ne=nei,nes)
      Rewind 1
      Write (6,190) (ecm(ne),(eta(nwave,ne),d(nwave,ne),nwave=1,nwaves),
     1ne=1,nes)
c     -------      calculate kinematics
      Do 60 ne=1,nes
         s = ecm(ne)*ecm(ne)
         sqarg=(s-(mn+mpi)*(mn+mpi))*(s-(mn-mpi)*(mn-mpi))/4./s
	 If (sqarg.le.0) sqarg=1.1
         kap(ne) = Sqrt(sqarg)
c     -----------compute alphas and convert deltas
         kk = kap(ne)
         epikap = Sqrt(mpi*mpi+kk*kk)
         enkap = Sqrt(mn*mn+kk*kk)
         ftot = -(epikap+enkap)/4./pi/pi/epikap/enkap
         factor = ftot/kk
         Do 50 nwave=1,14
            etadum = eta(nwave,ne)
            ddum = d(nwave,ne)*pi/180.
            If (etadum.le.0.) etadum = 1.
            rewave(nwave,ne) = 0.
            imwave(nwave,ne) = 0.
            If (nwave.gt.nwaves) GoTo 50
            rewave(nwave,ne) = .5*etadum*sin(2.*ddum)*factor
            imwave(nwave,ne) = .5*(1.-etadum*cos(2.*ddum))*factor
 50      Continue
 60   Continue
c     store unfoldef amplituded on unit 3,can fold latter/
      Rewind 3
      Write (3,205) nes
      Write (3,200) (kap(ne),(rewave(nwave,ne),imwave(nwave,ne),nwave=1,
     114),ne=1,nes)
      Rewind 3
      ndata = 0
      Write (6,210) nes,nwaves
      Write (6,230) (ecm(ne),kap(ne),(rewave(nwave,ne),imwave(nwave,ne),
     1nwave=1,14),ne=1,nes)
      Do 90 ne=1,nes
         Do 70 nwave=1,14
            t(nwave) = rewave(nwave,ne)
            tx(nwave) = imwave(nwave,ne)
 70      Continue
c     form partial wave decomp of pi+p,pi+n amplit udes
         Call mixup (bn,bnf,t,tx,nifty(1))
c     temp storage for printout
         Do 80 nwave=1,14
            d(nwave,ne) = bn(nwave)
            If (nwave.le.12) eta(nwave+2,ne) = bnf(nwave+4)
 80      Continue
         eta(1,ne) = bn(15)
         eta(2,ne) = bn(16)
 90   Continue
      rhl1 = 0.
      Write (6,240)
      Write (6,250) (kap(ne),(d(nwave,ne),nwave=1,13,2),eta(1,ne),(d
     1(nwav,ne) ,nwav=2,14,2) ,eta(2,ne),ne=1,nes)
      Write (6,260)
      Write (6,250) (kap(ne),rhl1,rhl1,(eta(nwave,ne),nwave=3,13,2),rhl1
     1,rhl1,(eta(nwav,ne),nwav=4,14,2),ne=1,nes)
c     -----amplitudes all stored, interpolate on kappa,using 2 point lag
c     amplitudes all stored on 3, Read back
 100  Continue
      ndata = 0
      Read (3,205) nstore
      Read (3,200) (kap(ne),(rewave(nwave,ne),imwave(nwave,ne),nwave=1,
     114),ne=1,nstore)
      Rewind 3
      If (nes.le.nstore) GoTo 110
      Write (6,220) nstore
      nes = nstore
 110  Continue
      x = kappa
      If (kappa.gt.kap(nes)) x = kap(nes)
      If (kappa.gt.kap(1)) GoTo 130
c     use scattering lengths to calc low e delts, only s+p ne 0,eta==1
      Do 120 nwave=1,14
         Call tpirsl (kappa,nwave,ddum,w,tr,ti)
         t(nwave) = tr
         tx(nwave) = ti
 120  Continue
      GoTo 140
c     here t and tx refer to interpolated values of rewave and imwave
 130  Call lagrng (x,kap,t,rewave,nes,14,2,100,14)
      Call lagrng (x,kap,tx,imwave,nes,14,2,100,14)
 140  Continue
c     check for set of 2 body amplitudes to zero
      If ((nifty(7).ne.1).and.(nifty(7).ne.2)) GoTo 160
c     nifty(7)=1, set p, d, f.. waves to 0   n(7)=2, only d,f =0
      ni = 5
      If (nifty(7).eq.2) ni = 9
      Do 150 nl=ni,16
         bn(nl) = 0.
         bnf(nl) = 0.
 150  Continue
      Return
 160  Call mixup (bn,bnf,t,tx,nifty(1))
      Return
c     -----formats
c
 170  format (13f6.1/8f6.1/8f6.1)
 180  format (49h0 input pion-nucleon phase shifts ecm,(eta,delta),  
     1  19h see tpicm for code,/,34h If folding,25 e s inserted before,
     2  50h first Read in value with rowe,salomon,Landau fits ,/,
     3  12h nes nwaves=,2i5,//)
 190  format (1h ,13f8.2/8f8.2/8f8.2)
 200  format (7e11.4/6e11.4/6e11.4/6e11.4/4e11.4)
 205  format(i5)
 210  format (55h non folded eigenchannel tmatrix as stored on tape3 for
     1,i5,8henergies/26h rel and im parts for all ,i3,9h channels/)
 220  format (25h ********nes .gt. nstore=,i5,21hnes set =nstore *****,/
     1/)
 230  format (9h,ecm kap=,2f10.2/14x,8e11.4/14x,10e11.4/14x,10e11.4)
 240  format (103h0separate pi**plus*  proton  and neutron partial wave
     1---non-flip ---( t matrices at c.m. momemtum kap ,/41h 1st line re
     2al   parts, second line imag /80h kap   swave p    swaven    pwave
     3 p    pwave n dwave p dwave n fwafe p fwave n  //)
 250  format (f8.2,10x,8e12.4,/18x,8e12.4)
 260  format (104h0separate pi**plus*  proton  and neutron partial wave
     1--- spin-flip--- t matrices at c.m. momemtum  kap ,/40h 1st line r
     2eal  parts, second line imag /80h kap   swave p    swaven    pwave
     3 p    pwave n dwave p dwave n fwafe p fwave n  //)
      End
