      program lpott
c    tape1= pin phases,  tape2= g(p)'s  tape3= pin t's,  tape4=denom
c    tape5= run param    tape6=output   tape7=scratch    tape8=w.f. storage
c ** see r h landau, program lpott, 1981, modified for ridge/vax
c	helium  3 ff updated too
c      calculation  of  r-matrix  and  t-matrix  by  matrix inversion wi
c       potential in p space for spin o and 1/2, any energy
c              energy,momentum,mass in  mev units throughout prog
c     nifty(1)     =1,piminus
c                  =0,piplus
c                  =2, pi+,pi-(+fcoul)
c                    n.b. same strong if nift=2,different if=3
c                  =3, cex after pi+,pi-(+fcoul for +-)
c                  =4, pi 0 (fcoul turned off automatically)
c                  =5  , protons
c                  =6 ,cex with pi+ and pi0
c                  =7,    k+
c                  =8,    k0
c                  =9,    neutron
c                  =-1,    (k+,k0) cex after k+,k0
c                  =-2,    (p,n) cex after p,n
c     nifty(2)*5mev=binding e(-5.*n(2) used) for e3b
c     nifty(3)= 0, no w.f. in coord space
c               1. print out wf info on tape8
c               2: wf plus pi-nucleon t matrix
c               3:pi-nucleon tmatrix only
c               4: single scatt (when n(6)used else)
c               5:set t=r (nonunitarize) but with ms
c               6:double scattering for 1/2*1/2

c     nifty(4)= 0 sig=1, constant off-shell with on shell valus
c               1 nlsp form off shell
c               2 local laplacian form ,genrlzed to d $ f waves

c     nifty(5)= 0  use softh(k0)
c               1   sin(theat)
c               2   so(on shell),no angle tf
c               5  e3b energy
c               6  e3b energy with folded tmatrice
c               7   e3b with aay magic vector prescription for th
c               9  soth(ko), aay

c     nifty(6)=  0,  nospin flip and no calc of ds and ts (0 x 0)
c                1  % calculate r matrix up to and include ds
c                2 % calculate r matrix up to and include ts
c                3  % include spin flip terms (0 x 1/2)
c                4  no scatt (u=0) as check of wf
c                5  r matrix for only ss
c     nifty(7)  =0    no true absp, all 2-body partial waves
c               =1     no true absp, only s(l=0) pi-n waves
c               =2     no true absp,  only s&p(l=o,1) pi-n waves
c               =3   true absp (old)
c               =4   true absp(new,yoo-landau model)

c     nifty(8)  =0, no fold,no read in phases, use stored t
c               =1,  fold + stored t matrix, dont read in phases
c               =2,  read in phases, dont fold
c               =3,  read in phases, calc t matrix,fold these t s
c               =4 for k+, use martins analytic form
c               =5,  k+,martin s analytic, fold

c     nifty(9)  =0   reduced output
c               =1   full output
c               =2 reduced output,but including pauli

c     nifty(10) =0   no coulomb amplitude
c                1  coulomb amplitude and phase change included
c               =2 coulomb amplitude, but not phase change
c               =3 exact coulomb via matching in momentum space
c               =4 exact coul + realistic form factor
c     nifty(11)= upwards shift in pi-(or k+) subenergy(pi-nucleon)
c          (12)                   pi+ (k0)
c          (13)                   pi-(or k+) nucleus channel energy
c          (14)                   pi+(k0)  nucleus  channel energy
c     nifty(15) =0     divide out proton ff from he ff
c                1     dont divide out proton ff

c     nifty(16) =0: modified gaussian rho and rho**2 & for p* z=n  nucleus will
c                   also use t(spin)*rho(matter)
c               =1: wood saxon rho,..but gaussian rho**2
c               =2: wood saxon rho and rho**2
c               =3: special form for he or d (ignored if n(19) ne 0)
c               =4: he3..fnsp=fnmat,fpsp=0
c               =5: he3..fnsp=fnmat=fpmt, fpsp=0
c               =6: use fmatter*tspin for 0*.5 (non-gaussian nucleus)
c               =7: uspin=(tsp*fsp+tspn*fmt) for 0*.5

c      nifty(17)=0: kmt potential
c               =1: watson potential,no a/a-1
c     nifty(18) =0: dont turn off any 1/2*1/2 amps
c               =1: no e terms
c               =2: no a+b terms
c               =3: no e and no a+b terms
c      nifty(19)=0: nifty(16) controls ff's
c               =1: hadjimichael form factors - no meson exchange currents
c               =2: hadjimichael form factors - with meson exchange currents
c      nifty(20)=0: no output to tape10 (for010.dat for vax)
c               =1: theta, ds/dw, polar, and r (no r for he3) to tape10

c *** nsex =0-normal/ =1-1st time thru/ =2-2nd time thru/ =3 diff cex

c *** plot up lab cross sections if nang was read in as a negative numb

c *** u(nk,1)= uplus for spin flip,  u(nk, 2)= uminus
c *** set = for nonflip (same code for tr,ti,trc,tic,tbr,tbi)
c       ----------------------------------------------
      implicit real*8  (a-h,o-z)
      double complex h,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,
     1           x13,x14,x15,x16,x17,x86,x76,x18,x19,x20
      complex*16 zi,ztan
      real*8 kapg,mn,mn2,k0,ko,ko2,k,k2,kk,mpi,mpi2,mr,kone
      real*8 ktwo,imtij,bnuc,bn,retij,xbarc,xi,xpi,xn,xgam,zcm,zch,theta
     1dsig,bnucf,bnf,ff,dflip,dnof,wsp,wsn,achp,achn,acmp,acmn,ymin1,
     2ymin2
c      nb need match dimen and no entries in data stm for ibm
      dimension m1(100), m2(100), y(50), gin(14,100), h(12,20), kapg(100
     1), tbr(25,4), tbi(25,4), tr(25,4), ti(25,4),x1(12),x2(1),x3(4),
     2 x4(3),x5(10),x6(7),x7(6),x8(6),x9(2),x10(5),x11(1),
     3 x12(1),x13(1),x14(1),x15(1),x16(4),x17(2),x18(4),x19(3),x20(2)
c      dimension of u/f depends on ngp (only need change in main)
c       double size for spin 0 case,so can up n grid pts
c      u(n,2) ,n=4*(ngp+1)**2 ,only need card needs cahanging
c      16pts-1156,24pts-2500, 32pts-4356, with other  obvious ones
c      nb also change equivalence of chart,as u and f enlarge
      dimension kk(66), u(4356,2), den(66), f(4356), gp(32), wt(32)
      dimension theta(95), yy(50), ff(4), dflip(95), dnof(95)
      dimension nifty(20), bnuc(20), bn(16), retij(14), imtij(14), bnucf
     1(20), dsig(95), chart(61,61), polar(61), denom(14,25), ecmden(25),
     2 bnf(16), fc(30), gc(30), fpc(30), gpc(30), scoul(30), fc0(30),
     3 gc0(30), fpc0(30), gpc0(30), reul(25,2), aimul(25,2)
      common /nlspfl/ gin,kapg,denom,ecmden,fmn,rew(14,100),aimw(14,100)
     1,akapd(100),nkapgs,necms,dummy(3934)
      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
      common /absp/ b0r,b0i,c0r,c0i
	common /nocheck/ nocheck
D      equivalence (u(1),y(1)), (u(51),yy(1)), (m2(1),dsig(1)), (chart(1,
D     11),u(1,1)), (gin(1,1),u(1,1)), (aimw(12,34),dnof(1)), (aimw(14,41)
D     2,dflip(1)), (aimw(2,49),polar(1)), (gin(1),reul(1)), (gin(51),
D     3aimul(1)),(x1(1),h(1,1)),(x2(1),h(1,2)),(x3(1),h(1,3)),
D     4(x4(1),h(1,4)),(x5(1),h(1,5)),(x6(1),h(1,6)),(x7(1),h(1,7)),
D     5(x8(1),h(1,8)),(x9(1),h(1,9)),(x10(1),h(1,10)),(x11(1),h(1,11)),
D     6(x12(1),h(1,12)),(x13(1),h(1,13)),(x14(1),h(1,14)),(x15(1),h(1,15)
D     7),(x16(1),h(1,16)),(x17(1),h(1,17)),(x18(1),h(1,18)),
D     8(x19(1),h(1,19)),(x20(1),h(1,20))
c      for ibm need remove (u,h) equivalence
	equivalence 
     3(x1(1),h(1,1)),(x2(1),h(1,2)),(x3(1),h(1,3)),
     4(x4(1),h(1,4)),(x5(1),h(1,5)),(x6(1),h(1,6)),(x7(1),h(1,7)),
     5(x8(1),h(1,8)),(x9(1),h(1,9)),(x10(1),h(1,10)),(x11(1),h(1,11)),
     6(x12(1),h(1,12)),(x13(1),h(1,13)),(x14(1),h(1,14)),(x15(1),h(1,15)
     7),(x16(1),h(1,16)),(x17(1),h(1,17)),(x18(1),h(1,18)),
     8(x19(1),h(1,19)),(x20(1),h(1,20))
c      for ibm need remove (u,h) equivalence
c
c
      data hbarc,mpi,pi/197.3289e+0,139.578e+0,3.141593e+0/,x1    /
     1 10h(p,n)     ,10h(k+,ko)   ,10hpi+       ,10hpi-       ,
     2 10hpi+,pi-   ,10hpi cex    ,10hpio       ,10hproton    ,
     3 10hcex(+pio) ,10hk+        ,10hk0        ,10hneutron   /,
     4x2     /10h*5=be         /,x3   /10hno w.f.    ,10hwave funct,
     5 10hw.f.+tpin  ,10htpi       /,x4    /10htpin=const,
     6 10hnlsp-g(p)  ,10hlocal lapp  /,x5    /10hso,kapo    ,
     7  10hsin,kap0   ,10hso,no ang  ,2*10hundef     ,10he3b,no aay
     8,10hundef      ,10he3b, aay   ,10hundefined ,10hso,aay    /,x6/
     910hno spin ms,10hss+ds      ,10hss+ds+ts   ,10hspin,ms    ,
     a 10hus=uc=0   ,10hundef     ,10hus=0,ucne0 /,x7    /
     b 10hno abs/s-f,10hno abs/s-s,10hno abs/s-p  ,10hvabs-old  ,
     c10hvabs-yool  ,10habs-old,sp/,x8    /10hno del/fld,10hfld,no del ,
     d10hdel,no fld,10hdel in,fld,10hk+ martin ,10hmartin,fld /,x9/
     e 10hno pauli  ,10hpauli tpin /,x10/10hno coul    ,10h+fc+phase ,
     f10h+fc        ,10hexact coul ,10hexact+f(q) /,x11    /
     g10hpi0,pi-,k+/,x12    /10hpi+,ko     /,x13    /
     h 10hpi-channel  /,x14    /10hpi+channel/,x15    /
     i10hnot used   /,x16    /10hrhog,rho2g,10hrhows,rh2g,10hrhows,r2ws,
     j10hhe/d      /,x17    /10hkmt       ,10hwatson    /
      data x18/10hfull amps ,10hno e      ,10hno a+b    ,10hno e,a,b   /,
     1 x19/10huse nif-16,10hhadjm's ff,
     2 10hhadj + mec/
     3 ,x20/10hno t10 plt,10htpe10 plot/
c
      areal(xx,yx) = -xx/(1+xx*xx+yx*yx-2*yx)/((ko/hbarc)**(2*ldum-1))
      aim(xx,yx) = (xx*xx+yx*yx-yx)/(1+xx*xx+yx*yx-2*yx)/((ko/hbarc)**(2
     1*ldum-1))
c
	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')
ccall ctchsig()
	write(6,879)
 879	format(' pion scattering via lpott(vaxversion) on 32,1/86')
cstop 'check'
	write(*,*) 'ok tape6'
      ztan = (0.e+0,0.e+0)
c     cp0=cpu_time (dumm)
 10   continue
c     mn here refers to nucleus mass,xn (in common) to nucleon
c     read  momenta  in  centre _of_mass
c      nr.le.0  is  nonrelativistic  nr .gt. 0 >     relativistic  case
c                                    .gt. 4 for approximate klein-gordon
      read (5,*,end=650) nr,lxmax,nocheck
      write (6,740) nr,lxmax
c     -------------    read  in  lab  kinetic  energy
      read (5,*) e
      write (6,770) e
      read (5,*) ngp,kode,b,nang,ymin1,ymin2
      if (nang.eq.0) nang = 3
      if ((ymin1+ymin2).ne.0.) go to 20
      ymin1 = -2.
      ymin2 = -4.
 20   continue
      nplot = 0
      if (nang.lt.0) nplot = 1
      nang = iabs(nang)
      write (6,750) ngp,kode,b,nang,ymin1,ymin2
      if (ngp.le.0) go to 650
      write (6,790) ngp
c     -----read in data
      read (5,*) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut
      write (6,880) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut
      read (5,*) nz,na,(nifty(n),n=1,20)
      write (6,830) nz,na,(nifty(n),n=1,20)
c     print out nifty
      do 25 nif=1,20
      jj=0
         j = nifty(nif)+1
         if (nif.eq.1) j = nifty(nif)+3
         if (nif.eq.2.or.nif.eq.11.or.nif.eq.12.or.nif.eq.13.or.nif.eq.
     1   14 .or.nif.eq.15)jj=1
         if (jj.eq.1) j=1
         write (6,710) nif,nifty(nif),h(j,nif)
         if (jj.eq.1) write(6,711)
 25   continue
      read (5,*) nes,nwaves,b0r,b0i,c0r,c0i
      write (6,910) nes,nwaves,b0r,b0i,c0r,c0i
      if (b0i.gt.1.e-03.or.c0i.gt.1.e-03) go to 30
c      use old values if nothing is read in
      b0r = -0.04
      b0i = 0.04
      c0r = 0.
      c0i = 0.08
      write (6,910) nes,nwaves,b0r,b0i,c0r,c0i
 30   if (na.ne.3) go to 40
c      square input size params for he3(prevoius read in squaed)
      achp = achp*achp
      acmp = acmp*acmp
      wsp = wsp*wsp
c	mod for neutron sizes too
	achn=achn*achn
	acmn=acmn*acmn
	wsn=wsn*wsn
      write(6,880)achp,acmp,wsp,achn,acmn,wsn
 40   continue
c       added for wf
      if (nifty(3).eq.1) rewind 8
      if ((nifty(1).eq.7).or.(nifty(1).eq.8).or.(nifty(1).eq.-1)) mpi =
     1493.668
      amass = na
      zcm = acmp
      zch = achp
      zws = wsp
      nsex = 0
      nifty1 = nifty(1)
c     for pi 0(ko) , use pi0(ko) mass and no coul
      if (nifty1.eq.4.or.nifty1.eq.8) nifty(10) = 0
      if (nifty1.eq.4) go to 60
c     charge exchage, do both pi+ and pi-,then subtract
c     or pi+ and pi 0 ,or k+,ko
      n = nifty(1)+3
      go to (50,50,60,60,50,50,60,60,50,60,60,60), n
 50   nsex = 1
      nifty(1) = 0
      nifty(10) = 1
      if (nifty1.eq.-1) nifty(1) = 7
      if (nifty1.eq.-2) nifty(1) = 5
      if (nifty1.eq.2) write (6,670)
      if (nifty1.eq.3) write (6,660)
      if (nifty1.eq.6) write (6,680)
      if (nifty1.eq.-1) write (6,690)
      if (nifty1.eq.-2) write (6,700)
 60   continue
c         mn here refers to nucleus mass,xn is nucleon
      write(6,780) amass
      mn = 938.256*amass
      mn2 = mn*mn
      mpi2 = mpi*mpi
      xn = (nz*938.256+(na-nz)*939.55)/na
      xbarc = hbarc
      xi = pi
      xpi = mpi
      mr = mn*mpi/(mn+mpi)
c     shift channel energy by nifty13 for pi,nifty14 for pi+
 80   if (nifty(1).eq.1) e = e+nifty(13)-nifty(14)
      if (nifty(1).eq.4.or.nifty(1).eq.8.or.nifty(1).eq.9) e = e+nifty
     1(13)
      if (nifty(1).eq.0.or.nifty(1).eq.7.or.nifty(1).eq.5) e = e+nifty
     1(14)
      write (6,830)
      if (nr.le.0) write (6,1000)
c     now k0 for nr case is the same as for reltv, onlywave eqtn differs
      rhl1 = mn+mn
      el = e+mpi
      ecm = sqrt(mn2+mpi2+2.*mn*el)
      k0 = sqrt((el*el-mpi2)*mn2/ecm/ecm)
      pl = sqrt(el*el-mpi2)
      write (6,960) e,el,ecm,k0,pl
      write (6,980)
      if (nr.gt.4) write (6,990)
      write (6,940)
      write (6,950) e,k0
      n1 = ngp+1
      n2 = n1+n1
      a = k0
c       special check if b=0,use different distrb of points
      if (b.eq.0.) a = 200.
      call gauss (ngp,kode,a,b,gp,wt)
      piinv = 1./pi
      rhl1 = mr+mr
      rhl2 = (rhl1+rhl1)*piinv
      rhl3 = piinv+piinv
      rhl4 = rhl3+rhl3
      ko = k0
c     cp = cpu_time (dumm) -cp0
c     write (6,800) cp
      ko2 = ko*ko
      emnos = sqrt(ko2+mn2)
      epios = sqrt(ko2+mpi2)
      if (nr.eq.0) emnos = mn
      if (nr.eq.0) epios = mpi
      s = (emnos+epios)**2
      el = e+mpi
      pl = sqrt(el*el-mpi2)
      sum = 0.
      if (nr.gt.0) go to 100
c      -----------------------------------------------------------------
c     nonrelativistic      denominator
c      -----------------------------------------------------------------
      psfac = ko*rhl1
      do 90 i1=1,ngp
         k = gp(i1)
         kk(i1) = k
         kk(n1+i1) = k
         k2 = k*k
         den(i1) = k2*wt(i1)/(k2-ko2)*rhl2
         den(n1+i1) = den(i1)
         sum = sum+wt(i1)/(k2-ko2)
 90   continue
      kk(n1) = ko
      kk(n2) = ko
      den(n1) = -sum*rhl2*ko2
      den(n2) = -sum*rhl2*ko2
      go to 120
c        ---------------------------------------------------------------
c     relativistic    denominator
c        ---------------------------------------------------------------
 100  continue
      psfac = 2.*ko*epios*emnos/(epios+emnos)
      sum = 0.
      do 110 i1=1,ngp
         k = gp(i1)
         kk(i1) = k
         kk(n1+i1) = k
         k2 = k*k
         emn = sqrt(k2+mn2)
         epi = sqrt(k2+mpi2)
         den(i1) = rhl3*k2/(emn+epi-emnos-epios)*wt(i1)
         if (nr.gt.4) den(i1) = rhl3*k2/((k2-ko2)/(psfac/ko))*wt(i1)
         den(i1+n1) = den(i1)
         sum = sum+wt(i1)/(k2-ko2)
 110  continue
      kk(n1) = ko
      kk(n2) = ko
      den(n1) = -sum*rhl4*emnos*epios/(emnos+epios)*ko2
      den(n2) = den(n1)
 120  continue
      lclass = ko*1.33*(na**(.333333))/hbarc
      ldmax = lxmax
      if (lxmax.eq.0) ldmax = 2.5*lclass+2.5
      lborn = 2*lclass+1
      if (lborn.lt.4) lborn = 4
      if ((na.lt.12).and.(e.lt.250.)) lborn = lborn+1
      if ((na.gt.25).and.(e.lt.150.)) lborn = lborn+5
      if (lborn.gt.ldmax) lborn = ldmax
      xgam=0.
      if (nifty(10).eq.0) go to 160
c     include coulomb, set uop
      zp = 1
      nzp = nifty(1)+1
      go to (150,130,140,140,140,150,140,150,140,140), nzp
 130  zp = -1
      go to 150
 140  zp = 0.
 150  xgam = zp*nz/137.036/(ko*(emnos+epios)/emnos/epios)
      alph = (nz-2)/3.
      alph = alph/2./(2.+3.*alph)
      if (na.lt.4) alph = 0.
c     set up sigl
      rhoc = ko*rcut/hbarc
      lcmax = ldmax
      if (lcmax.eq.1) lcmax = 2
      eta = xgam
      call sigcl (eta,scoul,lcmax,30)
      sig0 = scoul(1)
      if ((nsex.eq.2).and.(nifty1.ne.3)) go to 430
 160  continue
c     read in separable potentials
      if (nifty(4).ne.1) go to 211
      read (2,850) nkapgs
      do 170 nk=1,nkapgs
           read (2,860,end=180) kapg(nk),(gin(nwave,nk),nwave=1,14)
c    .end of file check, this is machine dependent
 170  continue
      go to 190
 180  write (6,841) nk
      nkapgs = nk
 190  rewind 2
      rewind 4
      read (4,850) necms,fmn
       if(necms.le.0) go to 210
      do 200 nk=1,necms
         read (4,870,end=201)ecmden(nk),(denom(nw,nk),nw=1,14)
 200  continue
      go to 210
 201  write(6,841) nk
      necms=nk
 210  continue
      if (nifty(8).ne.2) go to 211
      write(6,851) nkapgs
      write(6,861) (kapg(nk),(gin(nwave,nk),nwave=1,14),nk=1,nkapgs)
      write(6,851) necms,fmn
      write(6,871) (ecmden(nk),(denom(nw,nk),nw=1,14),nk=1,necms)
 211  continue
c     do loop over l for this e and  state
      do 420 ldum=1,ldmax
         ld = ldum
         l = ldum-1
         write (6,1070) l
c        cp = cpu_time (dumm) -cp0
c        write (6,800) cp
         call optp (kk,n1,ldum,ldmax,lborn,u,f,theta)
c     do loop over two spin sates(if necessary)
c        cp = cpu_time (dumm) -cp0
c        write (6,800) cp
         do 410 nspin=1,2
            n = nspin
            if ((nspin.eq.2).and.(nifty(6).ne.3)) go to 380
c     test if l large enuf for born approx
            if (ldum.gt.lborn) go to 340
            do 240 i2=1,n2
               do 240 i3=1,n2
                  i = (i3-1)*n2+i2
                  f(i) = u(i,nspin)*den(i3)
c      chek if just want ds/triple scattering
                  if ((nifty(6).eq.1).or.(nifty(6).eq.2)) f(i) = -f(i)
c     special ss only calc
                  if (nifty(6).eq.5) f(i) = 0.
                  if (nifty(6).ne.2) go to 230
c     ts do loop over index k(present indices i and j)
                  term = 0.
                  do 220 ik=1,n2
                     iik = (ik-1)*n2+i2
                     ikj = (i3-1)*n2+ik
                     term = term+u(iik,nspin)*den(ik)*u(ikj,nspin)
 220              continue
                  f(i) = f(i)+term*den(i3)
 230              if (i2.ne.i3) go to 240
                  f(i) = f(i)+1.
 240        continue
c       minv  does  matrix-inversion.
            if ((nifty(6).eq.1).or.(nifty(6).eq.2)) go to 250
            call minv (f,n2,d,m1,m2)
            write (6,810) d
c           if (ldum.le.1) call time (3,3)
            go to 260
 250        nms = 1+nifty(6)
            write (6,820) nms
 260        rr = 0.
            ri = 0.0
            rr1 = 0.
            ri1 = 0.
c      calculation  of  r-matrix   (on  shell)
            irhl = n2*(n1-1)
            irhl1 = n2*(n2-1)
            do 270 j3=1,n2
               j = j3*n2
               i = j-n1
               i1 = j3+irhl
               j1 = j3+irhl1
               rr = rr+f(i)*u(i1,nspin)
               ri = ri+f(j)*u(i1,nspin)
               rr1 = rr1+f(j)*u(j1,nspin)
               ri1 = ri1+f(i)*u(j1,nspin)
 270        continue
            write (6,970) rr,ri,rr1,ri1
            rr = -psfac*rr
            ri = -psfac*ri
            write (6,930) rr,ri
c      calculation  of    t-matrix  (on - shell),
c      t normalized to exp(i del)*sin(del)
 280        cnst1 = ri*ri
            cnst2 = rr*rr
            cnst3 = 1.+cnst1-ri-ri
            cnst4 = cnst3+4.*ri
            aovera = na/(na-1.)
            if (nifty(17).eq.1) aovera = 1.
            if (na.eq.2) aovera = 1.
            tr(ld,n) = rr/(cnst4+cnst2)*aovera
            ti(ld,n) = (ri+cnst1+cnst2)/(cnst4+cnst2)*aovera
            write (6,1030) tr(ld,n),ti(ld,n)
            ar = areal(tr(ld,n),ti(ld,n))
            ai = aim(tr(ld,n),ti(ld,n))
            if (e.le.30.) write (6,1060) ar,ai
c      calculation  ofphase  shift+  eta
            eta = sqrt(1+4*(tr(ld,n)**2+ti(ld,n)**2-ti(ld,n)))
            delta = 0.5*atan2(2.*tr(ld,n),(1.-2.*ti(ld,n)))*180./pi
            if (eta.gt.1.005) write (6,1020)
            write (6,1010) eta,delta
            trc = tr(ld,n)
            tic = ti(ld,n)
            ar = areal(trc,tic)
            ai = aim(trc,tic)
            if (e.le.30.) write (6,1060) ar,ai
            if (nifty(10).eq.0) go to 300
            if (nifty(10).lt.3) go to 290
c     ---- --coulomb v included on poten
c     determine phase shift by matching
c
c     determine tan delta(sl) from calculated t*s
            zi = (0.e00,1.e00)
            ztan = (1-2*ti(ld,n))+zi*(2*tr(ld,n))
            njump = 1
            if (njump.eq.1) go to 350
 290        continue
            njump = 0
c     include coulomb phase factor in nuclear amplitude
            sigl = scoul(ldum)
            if (nifty(10).eq.2) sigl = 0
            rhl5 = cos(2*sigl)
            rhl6 = sin(2*sigl)
            tr(ld,n) = rhl5*trc-rhl6*tic
            ti(ld,n) = rhl6*trc+rhl5*tic
            write (6,1050) tr(ld,n),ti(ld,n)
 300        continue
c     output info for calc of wf in rspace
c     need tape8 declared for storage
            if ((nifty(3).ne.1).and.(nifty(3).ne.2)) go to 330
            if ((ldum.gt.1).or.(nspin.eq.2)) go to 310
c     first print out of wf info
            rewind 8
            write (8,920) nz,na,(nifty(n),n=1,20),e,xgam
            write (8,880) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut
            lstore = ldmax
            if (lborn.lt.ldmax) lstore = lborn
            write (8,920) ngp,lstore
            write (8,720) (kk(n),n=1,n1)
c     write out l,nspin,r and f-1 matrices
 310        write (8,920) ldum,nspin
            write (8,720) rr,ri,trc,tic
            do 320 j=1,n1
               jreal = (j-1)*n2+n1
               jim = jreal+n1
               jrchk = (j+n1-1)*n2+n1
               jichk = jrchk+n1
               write (8,720) f(jreal),f(jim),f(jrchk),f(jichk)
 320        continue
 330        continue
            if (ldum.gt.lborn) go to 390
c       born  approximation   ------------------
 340        irbr = (n1-1)*n2+n1
            irbi = n1*n2
            n = nspin
            rbr = -u(irbr,n)*psfac
            rbi = -u(irbi,n)*psfac
            tbr(ld,n) = rbr*aovera
            tbi(ld,n) = rbi*aovera
            etab = sqrt(1+4*(tbr(ld,n)**2+tbi(ld,n)**2-tbi(ld,n)))
            deltb = 0.5*atan2(2.*tbr(ld,n),(1.-2.*tbi(ld,n)))*180./pi
            write (6,1080)
            write (6,930) rbr,rbi
            write (6,1030) tbr(ld,n),tbi(ld,n)
            ar = areal(tbr(ld,n),tbi(ld,n))
            ai = aim(tbr(ld,n),tbi(ld,n))
            if (e.le.30.) write (6,1060) ar,ai
            write (6,1010) etab,deltb
            if (nifty(10).lt.3) go to 370
c      match for born
            ztan = (1-2*tbi(ld,n))+zi*(2*tbr(ld,n))
 350        ztan = zi*(1-ztan)/(1+ztan)
c     determine tandelta(nl) via match
c      reset coul functions
            if (ldum.ne.1) go to 360
            eta = 0.
            rhoc = ko*rcut/hbarc
            lcmax = ldmax
            if (lcmax.eq.1) lcmax = 2
            accur = 1.e-14
            step = 999
            nminl = 0
            call rcwfn (rhoc,eta,nminl,lcmax,fc0,fpc0,gc0,gpc0,accur,
     1      step)
c     correction of couloum for a/(a-1)factor in kmt prescription
            eta = xgam/aovera
            call rcwfn (rhoc,eta,nminl,lcmax,fc,fpc,gc,gpc,accur,step)
            call sigcl (eta,scoul,lcmax,30)
 360        wron1 = fc(ld)*gpc0(ld)-gc0(ld)*fpc(ld)
            wron2 = fc(ld)*fpc0(ld)-fc0(ld)*fpc(ld)
            wron3 = fc0(ld)*gpc(ld)-gc(ld)*fpc0(ld)
            wron4 = gc0(ld)*gpc(ld)-gc(ld)*gpc0(ld)
            ztan = (ztan*wron1+wron2)/(wron3+ztan*wron4)
            ztan = (1+zi*ztan)/(1-zi*ztan)
            trc = -zi*ztan/2
            tic = (1-ztan)/2
            eta = abs(ztan)
            delta = 0.5*atan2(2*trc,(1-2*tic))*180./pi
            write (6,1040)
            write (6,1030) trc,tic
            write (6,1010) eta,delta
            ar = areal(trc,tic)
            ai = aim(trc,tic)
            if (e.le.30.) write (6,1060) ar,ai
c     include coulomb phase factor in nuclear amplitude
            if (njump.eq.1) go to 290
            sigl = scoul(ldum)
            rhl5 = cos(2*sigl)
            rhl6 = sin(2*sigl)
            tbr(ld,n) = rhl5*trc-rhl6*tic
            tbi(ld,n) = rhl6*trc+rhl5*tic
 370        continue
            if (ldum.le.lborn) go to 390
c     use b.a. fo rhese l*s
            ri = rbi
            rr = rbr
            go to 280
c     no dependence, set tl+=tl-
 380        tr(ld,2) = tr(ld,1)
            ti(ld,2) = ti(ld,1)
            tbr(ld,2) = tbr(ld,1)
            tbi(ld,2) = tbi(ld,1)
c     charge exchange check
 390        if (nsex.eq.0) go to 410
            l = ldum
            np2 = n+2
            rhl6 = cos(2.*sigl)
            rhl7 = sin(2.*sigl)
            if (nsex.eq.2) go to 400
c     1st time thru, store t"s
            tbr(l,np2) = tbr(l,n)
            tbi(l,np2) = tbi(l,n)
            tr(l,np2) = tr(l,n)
            ti(l,np2) = ti(l,n)
            go to 410
c     2nd time thru, take difference and store for latter
 400        rhl5 = sqrt(2.)
            if (nifty1.eq.6.or.nifty1.eq.-1.or.nifty1.eq.-2) rhl5 = 1./
     1      rhl5
            tbr(l,np2) = (tbr(l,np2)-tbr(l,n))/rhl5
            tbi(l,np2) = (tbi(l,np2)-tbi(l,n))/rhl5
            tr(l,np2) = (tr(l,np2)-tr(l,n))/rhl5
            ti(l,np2) = (ti(l,np2)-ti(l,n))/rhl5
 410     continue
c    .end of do loof over spin
 420  continue
c    .end of do loop over partial waves(after remove coulomb)
 430  if(nsex.ne.3)go to 435
      write(6,740) nsex
      write(6,1030) (tr(ll,1),ti(ll,1),tr(ll,2),ti(ll,2),ll=1,ldmax)
 435  continue
      write (6,1090)
      if (nifty(6).eq.3) write (6,1110)
      cosold = 2.
      noangs = (180./nang)+1
c     calculate angular disrtbn
      do 520 iang=1,noangs
         theta(iang) = (iang-1)*nang
         x = cos(theta(iang)*pi/180.)
         call legpol (x,y,20)
         gbr = 0.e00
         gbi = 0.e00
         gr = 0.e00
         gi = 0.e00
         if (nifty(6).ne.3) go to 440
         call plprme (x,yy,20)
 440     continue
         sumbr = 0.
         sumbi = 0.
         sumr = 0.
         sumi = 0.
         do 450 ldum=1,ldmax
            rhl = ldum-1.
            sumbr = sumbr+(ldum*tbr(ldum,1)+rhl*tbr(ldum,2))*y(ldum)
            sumbi = sumbi+(ldum*tbi(ldum,1)+rhl*tbi(ldum,2))*y(ldum)
            sumr = sumr+(ldum*tr(ldum,1)+rhl*tr(ldum,2))*y(ldum)
            sumi = sumi+(ldum*ti(ldum,1)+rhl*ti(ldum,2))*y(ldum)
            if (nifty(6).ne.3) go to 450
            gbr = gbr+(tbr(ldum,1)-tbr(ldum,2))*yy(ldum)
            gbi = gbi+(tbi(ldum,1)-tbi(ldum,2))*yy(ldum)
            gr = gr+(tr(ldum,1)-tr(ldum,2))*yy(ldum)
            gi = gi+(ti(ldum,1)-ti(ldum,2))*yy(ldum)
 450     continue
         coscm = x
         tivt = -2.*(kk(n1)**2)*(1.-coscm)
         qtf = 0.
         if (tivt.ge.0.) go to 460
         qtf = sqrt(-tivt)
 460     continue
         rhl = (sumr*hbarc/kk(n1))
         rhl5 = (sumi*hbarc/kk(n1))
         rhl6 = 0.
         rhl7 = (hbarc/kk(n1))**2
         dflip(iang) = rhl7*(1.-x*x)*(gr**2+gi**2)
         dflipb = rhl7*(1.-x*x)*(gbr**2+gbi**2)
         if (na.eq.2) dflipb = dflipb*2./3.
         if (na.eq.2) dflip(iang) = dflip(iang)*2./3.
         dnof(iang) = rhl7*(sumr**2+sumi**2)
         dnofb = rhl7*(sumbr**2+sumbi**2)
         dsig(iang) = dflip(iang)+dnof(iang)
         dsigma = dsig(iang)
         dsigb = dflipb+dnofb
c     polarization
         if (dsig(iang).eq.0) go to 470
         polar(iang) = 2.*rhl7*sqrt(1.000001-x*x)*(sumi*gr-sumr*gi)/dsig
     1   (iang)
 470     if (nifty(10).eq.0) go to 510
c     include coulomb ampl + phase changes
         ff(1) = 1
         if (nifty(10).ge.3) go to 500
c      reset raddii in case equiv wiped them out
         acmp = zcm
         achp = zch
         wsp = zws
         if (na.eq.2) go to 480
         if ((na.le.16).and.(na.gt.4)) go to 490
         if (nifty(16).eq.0) go to 490
c     for ca or non-gaussian nucleus
         call ffact (-tivt,ff,1)
         go to 500
 480     qr = qtf*rcoul/hbarc
         if (qr.eq.0.) go to 500
         ff(1) = (3./qr/qr)*(sin(qr)/qr-cos(qr))
         go to 500
 490     ff(1) = (1.-alph*(qtf*zcm/hbarc)**2)*exp(-.25*(qtf*zch/hbarc)**
     1   2)
 500     if (theta(iang).eq.0) theta(iang) = 1.d-2
c    changed e-10 to d-2 to evade vax overflow problems
         phasec = sin(theta(iang)*pi/360.)
         fcoul = (-xgam*ff(1)/2./kk(n1)/phasec/phasec)*hbarc
         phasec = -2.*xgam*log(phasec)
         rhl6 = phasec+2*sig0
         phasec = phasec+2.*sig0
         rhl = (sumr*hbarc/kk(n1))+fcoul*cos(phasec)
         rhl5 = (sumi*hbarc/kk(n1))+fcoul*sin(phasec)
         dsigma = rhl*rhl+rhl5*rhl5+dflip(iang)
         rhl = (sumr*hbarc/kk(n1))
         rhl5 = (sumi*hbarc/kk(n1))
         sumbr = sumbr*hbarc/kk(n1)+fcoul*cos(rhl6)
         sumbi = sumbi*hbarc/kk(n1)+fcoul*sin(rhl6)
         rhl6 = fcoul/sqrt(dsig(iang))
         dsigb = sumbr*sumbr+sumbi*sumbi+dflipb
         dnof(iang) = rhl*rhl+rhl5*rhl5
         dnofb = sumbr**2+sumbi**2
c     calc of t.f. to lab from cm, see goldberg+watson
c
 510     ektwo = mn+pl*pl*mn*(1.-coscm)/s
         ktwo = sqrt(ektwo**2-mn2)
         ekone = el+mn-ektwo
         kone = sqrt(abs(ekone**2-mpi2))
         coslab = sqrt(abs(kone**2-ko2*(1.-coscm**2)))/kone
         if (coslab.gt.cosold) coslab = -coslab
         cosold = abs(coslab)
         dot = (mn2+mpi2+2.*ekone*ektwo-s)/2.
         cmtol = kone*kone*s/(pl*mn*(ektwo*kone-dot*ekone/kone))
c      convert cross sections from fm**2 to mb
         dsigma = dsigma*10.
         dsigb = dsigb*10.
         dsig(iang) = dsig(iang)*10.
         dnof(iang) = dnof(iang)*10.
         dnofb = dnofb*10.
         dflip(iang) = dflip(iang)*10.
         dflipb = dflipb*10.
         dsigl = dsigma*cmtol
         dsigbl = dsigb*cmtol
         thlab = acos(coslab)*180./pi
         write (6,1100) theta(iang),coscm,dsigma,dsigb,tivt,thlab,coslab
     1   ,dsigl,dsigbl,rhl,rhl5,dsig(iang)
         if (nifty(6).eq.3) write (6,1120) dnof(iang),dnofb,polar(iang),
     1   dflip(iang),dflipb
         dsig(iang) = dsigma
         if (nifty(6).eq.3) go to 520
c     set up for plot
         dnof(iang) = dsig(iang)
         dflip(iang) = dsigb
         dsig(iang) = dsigma
 520  continue
      sigel = 0.
      sigtot = 0.
      brnel = 0.
      brntot = 0.
c     total cross section summations
c       remove coulomb phase from  trhs
c     n.b. the t*s still contain the short range part of the coulmb inte
c       so they are not pure nuclear
      do 560 ldum=1,ldmax
         if ((nifty(10).eq.0).or.(nifty(10).eq.2)) go to 540
         sigl = scoul(ldum)
         rhl5 = cos(2*sigl)
         rhl6 = sin(2*sigl)
         do 530 nspin=1,2
            trc = rhl5*tr(ldum,nspin)+rhl6*ti(ldum,nspin)
            tic = rhl5*ti(ldum,nspin)-rhl6*tr(ldum,nspin)
            tr(ldum,nspin) = trc
            ti(ldum,nspin) = tic
 530     continue
 540     continue
         sigtot = sigtot+ldum*ti(ldum,1)+(ldum-1.)*ti(ldum,2)
         brntot = brntot+ldum*tbi(ldum,1)+(ldum-1.)*tbi(ldum,2)
         sigel = sigel+ldum*(tr(ldum,1)**2+ti(ldum,1)**2)+(ldum-1.)*(tr
     1   (ldum,2)**2+ti(ldum,2)**2)
         brnel = brnel+ldum*(tbr(ldum,1)**2+tbi(ldum,1)**2)+(ldum-1.)*
     1   (tbr(ldum,2)**2+tbi(ldum,2)**2)
         if (nsex.ne.2) go to 560
c     2nd pass on charge exch,redefine t
         l = ldum
         do 550 n=1,2
            np2 = n+2
            tbr(l,n) = tbr(l,np2)
            tbi(l,n) = tbi(l,np2)
            tr(l,n) = tr(l,np2)
            ti(l,n) = ti(l,np2)
 550     continue
 560  continue
      rhl = 10.*4.*pi*(hbarc/kk(n1))**2
      sigtot = rhl*sigtot
      sigel = sigel*rhl
      brnel = brnel*rhl
      brntot = brntot*rhl
      rhl = sigtot-sigel
      write (6,1130) sigtot,sigel,rhl,e,xgam,brntot,brnel
      if (nplot.ne.1) go to 580
c      plot up cross sections on semilog 5 cycle paper
      nspin = 0
      if (nifty(6).eq.3) nspin = 1
c     varriables for plotting titles
      chart(1,1) = e
      chart(2,1) = achp
      chart(3,1) = acmp
      chart(4,1) = wsp
      chart(5,1) = achn
      chart(6,1) = acmn
      chart(7,1) = wsn
      chart(8,1) = nz
      chart(9,1) = amass
      do 570 ii=1,20
         j = ii+9
         chart(j,1) = nifty(ii)
 570  continue
      call logplt (theta,dsig,dnof,noangs,nspin,chart,ymin1)
      if (nifty(6).ne.3) call logplt (theta,dnof,dflip,noangs,1,chart,
     1ymin2)
      if (nifty(6).ne.3) go to 580
      call logplt (theta,dsig,dflip,noangs,1,chart,ymin2)
      call logplt (theta,polar,dflip,noangs,-1,chart,ymin2)
 580  if (nsex.eq.0) go to 620
      if (nsex.eq.3) go to 620
      if (nsex.eq.2) go to 610
      nsex = 2
      nifty(1) = 1
      if (nifty1.ne.-1) go to 590
c     k+done, do ko
      nifty(1) = 8
      mpi = 497.67
      nifty(10) = 0
      go to 80
 590  if (nifty1.ne.-2) go to 600
c      proton done, now do neutron
      nifty(1) = 9
      mpi = 939.5731
      nifty(10) = 0
      go to 80
 600  continue
c     nsex = 1, one pass for cex done,now do other
      if (nifty1.ne.6) go to 80
c     pi+ done, now do pi0
      nifty(1) = 4
      mpi = 134.9626
      nifty(10) = 0
      go to 80
 610  nsex = 3
      if (nifty1.ne.3.and.nifty1.ne.-1.and.nifty1.ne.-2) go to 620
      nifty(10) = 0
      ymin1 = ymin1-1
      go to 430
 620  continue
c     wave function in rspace print out as test
      if ((nifty(3).eq.0).or.(nifty(3).eq.3)) go to 640
      nspwf = 1
      lmaxwf = 5
      rhl = 0.2
      do 630 nr=1,200
         rf = (nr-1.)*rhl
         if (rf.eq.0.) rf = 1.e-3
         rmev = rf/hbarc
         if (rf.gt.rcut*2.) go to 640
         call wavfn (nspwf,lmaxwf,rmev)
         write (6,890) rf,(reul(l,1),aimul(l,1),l=1,lmaxwf)
 630  continue
 640  continue
      go to 10
c     _______________     formats
 650  write(6,841)
      stop
c
 660  format (53h  pi-plus, pi-minus and then cex(pi- to pi0) 3 calcs )
 670  format (26h  pi-plus,piminus 2 clacs )
 680  format (54h   pi-plus, pi-zero(with different mass) and the cex  )
 690  format (38h k+ ,ko(with mass difference),then cex)
 700  format (15h p, n, then cex)
 710  format(7h nifty(,i2,2h)=,i3,10x,a8,a2)
 711  format(1h$,36x,5hshift)
 720  format (5e16.7)
 730  format (4x,4i4)
 740  format (1h ,5i4)
 750  format (1h ,i4,i3,e10.2,i4,2f5.2)
 760  format (i4,i3,e10.2,i4,2f5.2)
 770  format (f10.4)
 780  format (1h ,30x,17h********** amass=,f10.4,10h**********,//7x
     1 ,4h tpi,10x,6hepilab,8x,10hecm(sqrts),4x,3hkcm,11x,4hplab )
 790  format (1h ,18hno of grid points=,i4)
 800  format (1h0,12x,15h time(seconds)=,f10.2)
 810  format (1h ,13x,34hdeterminant from matrix inversion=,e14.7/)
 820  format (47h r matrix calculated with ms series up to order,i3)
 830  format (1h ,2i5,1x,i2,i2,8i1,10i3)
 840  format (1h ,i5,e12.4,10x,6e13.6)
 841  format(12h eof in main ,i3)
 850  format (i3,f15.5)
 851  format(29h g(p)-unit4/d(e)-unit2 print ,/i3,f15.5)
 860  format (f10.4/4(6f13.6,/),4f13.6)
 861  format(1x ,150(f10.4/6f13.6/6f13.6/2f13.6/))
 870  format (e10.4/8e10.3/6e10.3)
 871  format(1x ,150(f10.4/8e10.3/6e10.3/))
 880  format (8f10.4)
 890  format (1x,f10.3,10e10.3)
 900  format (2i5,4f10.5)
 910  format (28h nes,nwaves,b0r,b0i,c0r,c0i=,2i5,4f10.5)
 920  format (2i5,1x,i2,i2,8i1,10i3,f10.3,e13.6)
 930  format (1h ,26h  normalized r(k,k..k)=   ,11x,2e14.7)
 940  format (1h0,10x,8h  energy,10x,15h    momentum   ,/)
 950  format (1h ,7x,f14.4,10x,f14.4)
 960  format (1h ,5f14.4)
 970  format (1h ,37h  unnormalized r matrix(re,im,re,-im),4e14.7)
 980  format (1h0,31h      relativistic  calculation)
 990  format (1h ,10x,33happroximate klein gordon equation)
 1000 format (1h0,47h  nonrelativstc calc--but with same k0 as reltv,40x
     1,32h  i.e. only the wave eqtndiffers)
 1010 format (1h$,85x,4heta=,f10.4,8h  delta=,f10.4)
 1020 format (41h0********** eta greater than 1 **********,//)
 1030 format (1h ,15h  t(k,k,k)=    ,22x,4e14.7)
 1040 format (38h now with exact coulomb after matching)
 1050 format (20x,40hcoulomb phase modified nuclear amplitude/16h   t(k,
     1k.k) =   ,2e14.7)
 1060 format (1h ,21h scatt length(fermi)=,2e14.4)
 1070 format (1h0,30h----- angular   momentum  =   ,i4,5h-----)
 1080 format (1h ,8x,23h  born  approximation  )
 1090 format (19h0  th-c.m.  cos(cm),3x,10hdsig/dw-cm,2x,11hds/dw(b.a.),
     13x,1ht,4x,6hth-lab,6x,8hcos(lab),2x,5hds/dw,4x,11hds/dw(b.a.),2x,
     210href(th)-cm,1x,3himf,4x,6hsignuc,//)
 1100 format (1h ,f8.1,f10.3,3e11.3,f8.1,f10.3,2e11.3,2f9.3,e12.4)
 1110 format (1h+,21x,6hnoflip,7x,8hnoflip-b,5x,5hpolar,9x,4hflip,9x,
     16hflip-b,//)
 1120 format (1h ,18x,5e11.3)
 1130 format (1h0,22htotal cross-section  =,f10.3,9x,8helast , ,
     120hinel cross-section =,2f10.3//19h cross-sections  in,22h milliba
     2rns  ---------,25h energy(kinetic),gamcoul=,f10.1,e11.3/20h born a
     3pprox total =,f10.3,28h elastic              sigma=,f10.3)
      end
