      program mpfasesh
c     pi-nucl phase shifts
      Implicit Real*8 (a-h,o-z)
      Real*8 mn,mpi,imtij,xi,xpi,xn,bnuc,b,bnf,retij,ymin1,ymin2,
     1wsp,wsn,achp,achn,acmn      
      Dimension bnuc(20),bn(16),bnucf(20),bnf(16), retij(14),imtij(14),
     1 nifty(20)
      Common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,hbarc,pi,mpi,xn,nz,nes
     1,nwaves,nifty,na
      Common /sec4/ achp,acmp,wsp,achn,acmn,wsn/sec6/rcoul,rcut
      data hbarc,mpi,pi/197.3289e+0,139.578e+0,3.141593e+0/

      Open(1,file='tape1',status='unknown')
      Open(2,file='tape2',status='unknown')
      Open(3,file='tape3',status='unknown')
      Open(4,file='tape4',status='unknown')
      Open(7,status='scratch',access='sequential',
     1				form='unformatted')
      Open(6,file='fases.dat',status='unknown')
      Call readin(e)
      mn=na*938.256
      Write(6,*)'niftys ante ', (nifty(j), j=1,20)
      Write(6,*)'na  mn ', na,mn
      Call momentum(e,mn,mpi)
      Write(6,*)'niftys desp', (nifty(j), j=1,20)
      Close(6)
      Close(1)
      Close(2)
      Close(3)
      Close(4)
      Close(7)
      Stop
      End
      
      Subroutine readin(e)
c     Reads tape1, to initialize basic input for pi-He3     
c     Outputs to disk file "fases.dat"
c     Lines taken from lpott
      Implicit Real*8 (a-h,o-z)
      Real*8 mn,mpi,imtij,xi,xpi,xn,bnuc,b,bnf,retij,ymin1,ymin2,
     1wsp,wsn,achp,achn,acmn
      Dimension bnuc(20),bn(16),bnucf(20),bnf(16), retij(14),imtij(14),
     1 nifty(20)     
      Common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,xbarc,xi,xpi,xn,nz,nes
     1,nwaves,nifty,na
      Common /sec4/ achp,acmp,wsp,achn,acmn,wsn/sec6/rcoul,rcut
      data hbarc,mpi,pi/197.3289e+0,139.578e+0,3.141593e+0/
      Open(5,file='tape5',status='unknown')
c      Read (5,*,end=650) nr,lxmax
      Read (5,*) nr,lxmax
      Write (*,740) nr,lxmax
c     -------------    Read  in  lab  kinetic  energy
      Read (5,*) e
      Write (*,770) e
 770  format (f10.4) 
      Read (5,*) ngp,kode,b,nang,ymin1,ymin2
      If (nang.eq.0) nang = 3
      Write (6,*) ngp,kode,b,nang,ymin1,ymin2
c     -----Read in data
      Read (5,*) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut
      Write (*,750) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut
      Read (5,*) nz,na,(nifty(n),n=1,20)
      Write (*,830) nz,na,(nifty(n),n=1,20)
c     print out nifty
      Do 25 nif=1,20
          j = nifty(nif)+1
         If (nif.eq.1) j = nifty(nif)+3
ccc next line added with removal of nifty(8) = 0, 1
         If (nif.eq.8) j = nifty(nif)-1
   25   Continue
      Read (5,*) nes,nwaves,b0r,b0i,c0r,c0i
      Write (*,910) nes,nwaves,b0r,b0i,c0r,c0i
      Close(5)
c 650  Write(6,841)
 710  format(7h nifty(,i2,2h)=,i3,10x,a8,a2)
 740  format (1h ,5i4)
 750  format (1h ,i4,i3,e10.2,i4,2f5.2)
 830  format (1h ,2i5,1x,i2,i2,8i1,10i3)
c 841  format(12h eof in main ,i3)
 910  format (28h nes,nwaves,b0r,b0i,c0r,c0i=,2i5,4f10.5)
      Return
      End
      
      Subroutine momentum(e,mn,mpi)
c     computes kappa0 pi-nucleon com momentum  and calls tpicm
c           mn: nucleus mass
c           mpi
      Implicit Real*8 (a-h,i,k,m,o-z)
      mn2 = mn*mn
      mpi2 = mpi*mpi
      el = e+mpi
      ecm = Sqrt(mn2+mpi2+2.*mn*el)
      k0 = Sqrt((el*el-mpi2)*mn2/ecm/ecm)
      s = k0*k0
      If (s.ne.0.) kappa0 = Sqrt((s-mnpi2)*(s-mnmpi2)/4./s)
      Write(6,*)'momentum  energy ecm  ', kappa0, e, ecm
      Call tpicm(kappa0)
      Return
      End

      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 ampkitudes 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

      Subroutine tpirsl (kappa,ni,ddum,w,tr,ti)
c *** see RH Landau, Program lpott, 1981
c     calc of pin phases  and amplitudes with Rowe,Salomon,Landau
c     formula, or Landau s for higher l
      Implicit Real*8 (a-h,i,k,m,o-z)
      Dimension a(14), b(14), c(14), xx(14), w0(14), q0(14), gam(14)
c       use rowe,salomon,Landau for low e s,p,d13,d15;Landau for others
c
      data mpi,mn,pi,kold/139.578e+0,938.903,3.141593e+0,0./,
     1 a/-11.2,-2.91,11.4,16.8,-5.7,-1.31,11.2e-13,-3.15e-12,
     2 8.31e-18,5.46e-17,0.109,0.112,-38.26e-18,7.43e-18/,
     3 b/-30.7,3.45,-15.4,-35.4,25.4,1.22,-0.369e-14,0.118e-13,
     4 -0.363e-19,-0.236e-18,-0.031,-0.270,0.205e-18,-0.222e-19/,
     5 c/21.,-1.5,7.2,27.,-29.,-0.4,0.,0.,0.,0.,0.003,0.19,0.,0./,
     6 xx/.31,.22,.99,.44,.61,.23,0.,0.,0.,0.,.54,.43,0.,0./,
     7 w0/1655.,1850.,1233.,1550.,1435.,1815.,0.,0.,0.,0.,
     8 1525.,1670.,0.,0./,
     9 q0/550.,678.,228.,477.,393.,656.,0.,0.,0.,0.,459.,560.,0.,0./,
     a gam/170.,200.,116.,105.,230.,255.,0.,0.,0.,0.,125.,155.,0.,0./
      If (abs(kappa-kold).lt.1.e-03) GoTo 10
      kold = kappa
      rhl = kappa*kappa
      If (kappa.eq.0.) kappa = 1.e-10
      epikap = Sqrt(mpi*mpi+rhl)
      enkap = Sqrt(mn*mn+rhl)
      w = epikap+enkap
      s = w*w
      factor = (-1./kappa)*w/4./pi/pi/epikap/enkap
      q = kappa/mpi
      q2 = q*q
      q4 = q2*q2
c     determine l value for formula
 10   If ((ni.eq.1).or.(ni.eq.4)) l = 0
      If ((ni.eq.2).or.(ni.eq.3).or.(ni.eq.5).or.(ni.eq.6)) l = 1
      If ((ni.eq.7).or.(ni.eq.8).or.(ni.eq.11).or.(ni.eq.12)) l = 2
      If ((ni.eq.9).or.(ni.eq.10).or.(ni.eq.13).or.(ni.eq.14)) l = 3
      jj = 2*l+1
      If ((ni.eq.7).or.(ni.eq.8).or.(ni.eq.9).or.(ni.eq.10).or.(ni.eq.13
     1).or.(ni.eq.14)) GoTo 20
      rhl = q0(ni)/mpi
      del = (q**jj)*(a(ni)*1.e-2+b(ni)*q2*1.e-3+c(ni)*q4*1.e-4+xx(ni)*
     1gam(ni)*w0(ni)/(rhl**jj)/(w0(ni)*w0(ni)-s))
      ddum = atan(del)
      GoTo 30
c     no salomon fit, use old rhl fit
 20   ddum = (a(ni)+b(ni)*kappa)*(kappa**jj)*pi/180.
 30   tr = 0.5*factor*sin(2.*ddum)
      ti = 0.5*factor*(1.-cos(2.*ddum))
      Return
      End
      Subroutine mixup (bn,bnf,rt,ait,nif1)
c *** see RH Landau, Program lpott, 1981
c     mixup calls pion, nucleon or kaon mix to for partial waves
c     from eigen channel amplitudes
      Implicit  Real*8   (a-h,o-z)
      Dimension bn(16), bnf(16), rt(14), ait(14)
      n = nif1+1
      GoTo (10,10,10,10,10,30,10,20,20,30), n
 10   Call pimix (bn,bnf,rt,ait)
      Return
 20   Call kmixup (bn,bnf,rt,ait)
      Return
c     protons not in this version
 30   STOP
      End

      Subroutine pimix (bn,bnf,t,tx)
      Implicit  Real*8   (a-h,o-z)
c *** see RH Landau, Program lpott, 1981
      Dimension bn(16), bnf(16), t(14), tx(14)
c     form partial wave decomp of pi+p,pi+n amplit udes
c     -non-flip, sp,sp,pp,pn,...
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    7 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
      rhl = 2./3.
      bn(1) = t(1)
      bn(2) = tx(1)
      bn(3) = bn(1)/3.+rhl*t(4)
      bn(4) = bn(2)/3.+rhl*tx(4)
      bn(5) = t(2)+2.*t(3)
      bn(6) = tx(2)+2.*tx(3)
      bn(7) = bn(5)/3.+rhl*(t(5)+2.*t(6))
      bn(8) = bn(6)/3.+rhl*(tx(5)+2.*tx(6))
      bn(9) = 2.*t(7)+3.*t(8)
      bn(10) = 2.*tx(7)+3.*tx(8)
      bn(11) = bn(9)/3.+rhl*(2.*t(11)+3.*t(12))
      bn(12) = bn(10)/3.+rhl*(2.*tx(11)+3.*tx(12))
      bn(13) = 3.*t(9)+4.*t(10)
      bn(14) = 3.*tx(9)+4.*tx(10)
      bn(15) = bn(13)/3.+rhl*(3.*t(13)+4.*t(14))
      bn(16) = bn(14)/3.+rhl*(3.*tx(13)+4.*tx(14))
c     spon flip,same order as above
      bnf(1) = 0.
      bnf(2) = 0.
      bnf(3) = 0.
      bnf(4) = 0.
      bnf(5) = t(3)-t(2)
      bnf(6) = tx(3)-tx(2)
      bnf(7) = bnf(5)/3.+rhl*(t(6)-t(5))
      bnf(8) = bnf(6)/3.+rhl*(tx(6)-tx(5))
      bnf(9) = t(8)-t(7)
      bnf(10) = tx(8)-tx(7)
      bnf(11) = bnf(9)/3.+rhl*(t(12)-t(11))
      bnf(12) = bnf(10)/3.+rhl*(tx(12)-tx(11))
      bnf(13) = t(10)-t(9)
      bnf(14) = tx(10)-tx(9)
      bnf(15) = bnf(13)/3.+rhl*(t(14)-t(13))
      bnf(16) = bnf(14)/3.+rhl*(tx(14)-tx(13))
      Return
      End

      Subroutine lagrng (x,arg,y,val,ndim,nfs,nptx,maxarg,maxfs)
c *** see RH Landau, Program lpott, 1981
c     lagrange interpolation,unequally spaced points
c     npts=2,3,4,5,6, nfs functions(y*s) simultaneously interpltd
c     x= value of argument, arg is the tabulated x*s
c     y= a vector of interpolad functions, from tabulted val
c     ndim= Dimension of table
c     nfs= = of functions simult interpolated
c     maxarg and maxfs are maximum values of subscripts,ie Dimensions
      Implicit  Real*8   (a-h,o-z)
      Dimension arg(maxarg), val(maxfs,maxarg), y(maxfs)
c     -----find x0 the Closest point to x
      npts = nptx
c                                                                   2=st
      ni = 1
      nf = ndim
 10   If ((x.le.arg(ni)).or.(x.ge.arg(nf))) GoTo 30
      If ((nf-ni+1).eq.2) GoTo 70
      nmid = (nf+ni)/2
      If (x.gt.arg(nmid)) GoTo 20
      nf = nmid
      GoTo 10
 20   ni = nmid
      GoTo 10
c     ------ x is one of the tabltd values
 30   If (x.le.arg(ni)) GoTo 60
      nn = nf
 40   nused = 0
      Do 50 n=1,nfs
         y(n) = val(n,nn)
 50   Continue
      Return
 60   nn = ni
      GoTo 40
c     ------- 2 pts left, choose smaller one
 70   n0 = ni
      nn = npts-1
      GoTo (140,110,100,90,80), nn
 80   Continue
      If (((n0+3).gt.ndim).or.((n0-2).lt.1)) GoTo 90
      nused = 6
      GoTo 150
 90   Continue
      If ((n0+2).gt.ndim) GoTo 110
      If ((n0-2).lt.1) GoTo 100
      nused = 5
      GoTo 150
 100  Continue
      If (((n0+2).gt.ndim).or.((n0-1).lt.1)) GoTo 110
      nused = 4
      GoTo 150
 110  If ((n0+1).lt.ndim) GoTo 130
c     n0=ndim,special case
 120  nn = ndim
      GoTo 40
 130  nused = 3
      If ((n0-1).lt.1) nused = 2
      GoTo 150
 140  If ((n0+1).gt.ndim) GoTo 120
      nused = 2
 150  Continue
c     at least 2 pts left
      y0 = x-arg(n0)
      y1 = x-arg(n0+1)
      y01 = y1-y0
      c0 = y1/y01
      c1 = -y0/y01
      If (nused.eq.2) GoTo 160
c     at least 3 pts
      ym1 = x-arg(n0-1)
      y0m1 = ym1-y0
      ym11 = y1-ym1
      cm1 = -y0*y1/y0m1/ym11
      c0 = c0*ym1/y0m1
      c1 = -c1*ym1/ym11
      If (nused.eq.3) GoTo 180
c     ------at least 4 pts
      y2 = x-arg(n0+2)
      ym12 = y2-ym1
      y02 = y2-y0
      y12 = y2-y1
      cm1 = cm1*y2/ym12
      c0 = c0*y2/y02
      c1 = c1*y2/y12
      c2 = -ym1*y0*y1/ym12/y02/y12
      If (nused.eq.4) GoTo 200
c     at least 5 pts
      ym2 = x-arg(n0-2)
      ym2m1 = ym1-ym2
      ym20 = y0-ym2
      ym21 = y1-ym2
      ym22 = y2-ym2
      cm2 = ym1*y0*y1*y2/ym2m1/ym20/ym21/ym22
      cm1 = -cm1*ym2/ym2m1
      c0 = -c0*ym2/ym20
      c1 = -c1*ym2/ym21
      c2 = -c2*ym2/ym22
      If (nused.eq.5) GoTo 220
c     at least 6 pts
      y3 = x-arg(n0+3)
      ym23 = y3-ym2
      ym13 = y3-ym1
      y03 = y3-y0
      y13 = y3-y1
      y23 = y3-y2
      cm2 = cm2*y3/ym23
      cm1 = cm1*y3/ym13
      c0 = c0*y3/y03
      c1 = c1*y3/y13
      c2 = c2*y3/y23
      c3 = ym2*ym1*y0*y1*y2/ym23/ym13/y03/y13/y23
      GoTo 240
 160  Continue
      Do 170 n=1,nfs
         y(n) = c0*val(n,n0)+c1*val(n,n0+1)
 170  Continue
      GoTo 260
 180  Continue
      Do 190 n=1,nfs
         y(n) = cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0+1)
 190  Continue
      GoTo 260
 200  Continue
      Do 210 n=1,nfs
         y(n) = cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0+1)+c2*val(n,n0+
     1   2)
 210  Continue
      GoTo 260
 220  Continue
      Do 230 n=1,nfs
         y(n) = cm2*val(n,n0-2)+cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0
     1   +1)+c2*val(n,n0+2)
 230  Continue
      GoTo 260
 240  Continue
      Do 250 n=1,nfs
         y(n) = cm2*val(n,n0-2)+cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0
     1   +1)+c2*val(n,n0+2)+c3*val(n,n0+3)
 250  Continue
 260  Return
      End

      Subroutine kmixup (bn,bnf,rt,ait)
c *** see RH Landau, Program lpott, 1981
      Implicit  Real*8   (a-h,o-z)
      Dimension bn(16), bnf(16), t(28), rt(14), ait(14)
c
c     form partial wave decompositon of k+p, l+n ampli-
c     tudes. non flip sp,pp,pn ...etc p:proton,n: netron
c     code for bn,bnf
c     all for k plus
c     1 re s,p    2 im s,p    3 re s,n   4 im s,n
c     5 re p,p    6 im p,p    7 re p,n   8 im p,n
c     9 re d,p   10 im d,p   11 re d,n  12 im d,n
c     13re f,p  14 im f,p    15 re f,n   16 im f,n
c     spin non-flip
      Do 10 n=1,14
         nodd = 2*n-1
         neven = 2*n
         t(nodd) = rt(n)
         t(neven) = ait(n)
 10   Continue
      bn(1) = t(1)
      bn(2) = t(2)
      bn(3) = (t(1)+t(7))/2.
      bn(4) = (t(2)+t(8))/2.
      bn(5) = (2.*t(5))+t(3)
      bn(6) = 2.*t(6)+t(4)
      bn(7) = ((2.*t(5))+t(3)+(2.*t(11))+t(9))/2.
      bn(8) = ((2.*t(6))+t(4)+(2.*t(12))+t(10))/2.
      bn(9) = 3.*t(15)+2.*t(13)
      bn(10) = 3.*t(16)+2.*t(14)
      bn(11) = ((3.*t(15))+(2.*t(13))+(3.*t(21))+(2.*t(19)))/2.
      bn(12) = ((3.*t(16))+(2.*t(14))+(3.*t(22))+(2.*t(20)))/2.
      bn(13) = 3.*t(17)
      bn(14) = 3.*t(18)
      bn(15) = ((3.*t(17))+(3.*t(23)))/2.
      bn(16) = ((3.*t(18))+(3.*t(24)))/2.
c     spin flip same order as above
      bnf(1) = 0.
      bnf(2) = 0.
      bnf(3) = 0.
      bnf(4) = 0.
      bnf(5) = t(5)-t(3)
      bnf(6) = t(6)-t(4)
      bnf(7) = (t(5)-t(3)+t(11)-t(9))/2.
      bnf(8) = (t(6)-t(4)+t(12)-t(10))/2.
      bnf(9) = t(15)-t(13)
      bnf(10) = t(16)-t(14)
      bnf(11) = (t(15)-t(13)+t(21)-t(19))/2.
      bnf(12) = (t(16)-t(14)+t(22)-t(20))/2.
      bnf(13) = -t(17)
      bnf(14) = -t(18)
      bnf(15) = (-t(17)-t(23))/2.
      bnf(16) = (-t(18)-t(24))/2.
      Return
      End







