Subroutine wavfn (nspin,lmax,rmev)
c *** see RH Landau, Program lpott, 1981
c     wavfn calculates r space wf from r matrix and finverse
c     for all l s up to lmax at rmev (r in mev-1) and spin state nspin
c        ref(npt,l,nspin)  imf(npt,l,nspin)   the inverted matrix
c         rr(nspin,l)    ri(nspin,l)    the r matrix
c       reul(l,spin)  aimul(l,nspin)  the radial part of w.f.(withkr)
c     most arguments passed via Common nlspfl,with dummy filler        w
      Implicit Real*8(a-h,k,i,o-y)
      Implicit Complex*16(z)
      Common /nlspfl/ reul(25,2),aimul(25,2),ref(25,25,2),imf(25,25,2),
     1kk(34),rr(2,30),ri(2,30),fj(50),fc(30),gc(30),fpc(30),gpc(30),
     2znorm(25,2),tr(2,30),ti(2,30),nifty(20)
c
      data lcheck/-1/
      zi = (0.,1.)
      If (lcheck.gt.0) GoTo 100
c     first call, Read in all fmatrices and store them for latter r valu
      lcheck = 1
      Rewind 8
      Read (8,230) nz,na,(nifty(n),n=1,20),e,xgam
      Read (8,240) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut
      Read (8,230) ngp,lstore
      n1 = ngp+1
      If (lmax.gt.lstore) Write (6,210) lstore,lmax
      Read (8,250) (kk(n),n=1,n1)
      nspmax = 1
c     Do loop over l
      Do 40 l=1,lstore
         If (nifty(6).eq.3) nspmax = 2
c     Do loop over spin
         Do 30 nsp=1,nspmax
            znorm(l,nsp) = 0.
            Read (8,230) ldum,nspx
            Read (8,250) rr(nspx,l),ri(nspx,l),tr(nspx,l),ti(nspx,l)
            If ((nspx.ne.nsp).or.(ldum.ne.l)) Write (6,220)
c     Do loop over grid points
            Do 20 npt=1,n1
               Read(8,250,end=10)ref(npt,l,nsp),imf(npt,l,nsp),rhl1,rhl2
               If ((ref(npt,l,nsp).eq.rhl2).and.(imf(npt,l,nsp).eq.(-
     1         rhl1))) GoTo 20
 10            Write (6,220) nsp,l,npt,ref(npt,l,nsp),rhl1,imf(npt,l,nsp
     1         ),rhl2
               STOP
 20         Continue
 30      Continue
 40   Continue
c     calc wf normalization by matching to outgoing coul or spbesl
c     first calcl denom then  numerat
      Do 70 npt=1,n1
         rhl = kk(npt)*rcut/197.3286
         Call spbesl (lmax,rhl,fj)
         Do 60 l=1,lstore
            Do 50 nsp=1,nspmax
               zomega = na*(ref(npt,l,nsp)+zi*imf(npt,l,nsp))/
     1                     (1.-zi*(rr(nsp,l)+zi*ri(nsp,l)))/(na-1.)
               If (npt.eq.n1) zomega = zomega-1./(na-1.)
               znorm(l,nsp) = znorm(l,nsp)+fj(l)*zomega
 50         Continue
 60      Continue
 70   Continue
c     calc numerator
      eta = xgam
      If (nifty(10).ne.3) eta = 0.
      rhl = kk(n1)*rcut/197.3286
      lcmax = lmax
      accur = 1.e-14
      step = 999
      nmin = 0
      Call rcwfn (rhl,eta,nmin,lcmax,fc,fpc,gc,gpc,accur,step)
      Do 90 l=1,lstore
         Do 80 nsp=1,nspmax
            zidel = Sqrt(1-2*ti(nsp,l)+2*zi*tr(nsp,l))
            rhl1 = zidel
            rhl2 = zi*zidel
            zminus = rhl1+zi*rhl2
            zcos = (zidel+zminus)/2.
            zsin = (zidel-zminus)/2./zi
            znorm(l,nsp) = zidel*(zcos*fc(l)+zsin*gc(l))/rhl/znorm(l,nsp
     1      )
 80      Continue
 90   Continue
c     calculate wave function for l up to lstore,use jl l gt lstore
 100  If (rmev.ge.rcut/197.3286) GoTo 160
c     Do loop over grid points
      Do 150 npt=1,n1
         lCall = lmax+6
         Call spbesl (lcall,rmev*kk(npt),fj)
c     Do loop over spin and partial waves
         Do 140 nsp=1,nspmax
            Do 130 l=1,lmax
               If (l.gt.lstore) GoTo 120
               If (npt.ne.1) GoTo 110
c     intialize
               reul(l,nsp) = 0.
               aimul(l,nsp) = 0.
 110           zomega = na*(ref(npt,l,nsp)+zi*imf(npt,l,nsp))/
     1           (1.-zi*(rr(nsp,l)+zi*ri(nsp,l)))/(na-1.)
               If (npt.eq.n1) zomega = zomega-1./(na-1.)
               If (nifty(3).eq.2) zomega = zomega*(na-1.)/na
               If ((nifty(3).eq.2.).and.(npt.eq.n1)) zomega = zomega+1./
     1         na
               zomega = zomega*znorm(l,nsp)*fj(l)
               reul(l,nsp) = reul(l,nsp)+zomega
               aimul(l,nsp) = aimul(l,nsp)-zi*zomega
               GoTo 130
c      use undistored wabes for l gt lstore
 120           If (npt.ne.n1) GoTo 130
               reul(l,nsp) = fj(l)
               aimul(l,nsp) = 0.
               If (eta.eq.0.) GoTo 130
c     use Coulomb waves
               rhl = kk(n1)*rmev
               lcmax = lmax
               Call rcwfn (rhl,eta,nmin,lcmax,fc,fpc,gc,gpc,accur,step)
               reul(l,nsp) = fc(l)
 130        Continue
 140     Continue
 150  Continue
c     formats
      Return
c     use asympt form of waves for large r(gt rcut)
 160  rhl = kk(n1)*rmev
      lcmax = lmax
      Call rcwfn (rhl,eta,nmin,lcmax,fc,fpc,gc,gpc,accur,step)
      Do 200 l=1,lmax
         Do 190 nsp=1,nspmax
            If (l.le.lstore) GoTo 170
c      l gt lsore, use undistored waves
            zcos = 1.
            zsin = 0.
            zidel = 1.
            GoTo 180
 170        zidel = Sqrt(1-2*ti(nsp,l)+2*zi*tr(nsp,l))
            rhl1 = zidel
            rhl2 = zi*zidel
            zminus = rhl1+zi*rhl2
            zcos = (zidel+zminus)/2.
            zsin = (zidel-zminus)/2./zi
 180        zminus = zidel*(zcos*fc(l)+zsin*gc(l))/rhl
            reul(l,nsp) = zminus
            aimul(l,nsp) = -zi*zminus
 190     Continue
 200  Continue
      Return
c
 210  format (55h in wavfn lmax gt lstore,use plane wave for l gt lstore
     1,2i5)
 220  format (37h -----file Read error in wavfn or eof/39h       nsp,l,n
     1pt,refin,rhl1,imfin,rhl2=/3i5,4e20.7)
 230  format (2i5,2x,i1,i2,8i1,10i3,f10.3,e13.6)
 240  format (8f10.4)
 250  format (5e16.7)
      End