program psipolar
c     computes psi(r) for different values of the angle
c     theta see eq. 150 in lpott\  
c     this run is for CARBON 12 
c     for other nucleus change rcoul and rcut
      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)
      Dimension y(50),zu(20)
      Open(8,file='tape8',status='unknown')
      Open(9,file='psipolar.dat',status='unknown')
      pi=3.141592
      hbarc=197.329
      rcoul=3.0  !for carbon
      rcut=7.0    !for carbon
      nspin=1 
      lmax=5
      h=0.4 !step in r
      zi=(0.0, 1.0)
      theta=0.0 !degrees
      n=lmax+2
      Do ny=1,33              ! for y axis
          yy=(ny-17)*h
        Do nx=1,33            ! loop for x axis
          xx=(nx-17)*h
          radius=sqrt(xx**2+yy**2)
          if(xx.eq.0.0)theta=0.0
          if(xx.ne.0.0)theta=atan2(yy,xx)
          rad=cos(theta)
          Call legpol(rad,y,n)
          If(radius.eq.0.0)radius = 1.e-3
          rmev = radius/hbarc
c          If (radius.gt.rcut*2.) GoTo 10
          Call wavfn (nspin,lmax,rmev)
          Do j=1,lmax 
            zu(j)=(reul(j,1),aimul(j,1))
          End Do
          zsuml=(0.0,0.0)
          Do j=1,lmax           ! loop for L
            zsuml=zsuml+zu(j)*(zi**(j-1))*(2.0*(j-1)+1)*y(j)
          End Do
          psi2=real((zsuml))**2+(aimag(zsuml))**2
          Write(9,100)xx,yy,psi2
        End Do                    ! x
 10     Continue
        Write(9,*)'***'
        End Do                   ! y
 100  Format(1h ,3e13.5)
      Close(8)
      Close(9)
      Stop
      End

      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
           write(6,*)'ngp lstore',ngp, lstore
      n1 = ngp+1
      If (lmax.gt.lstore) Write (6,210) lstore,lmax
      Read (8,250) (kk(n),n=1,n1)
          write(6,410)(kk(n),n=1,n1)
 410       format(1h ,5e16.7)
      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
             write(6,*)'ldum,nspx',ldum,nspx   !mpaez
             write(6,*)'nspx nsp ldum l',nspx,nsp,ldum,l
            Read (8,250) rr(nspx,l),ri(nspx,l),tr(nspx,l),ti(nspx,l)
        write(6,400) rr(nspx,l),ri(nspx,l),tr(nspx,l),ti(nspx,l)
 400    format(1h , 4(e16.7))
             If ((nspx.ne.nsp).or.(ldum.ne.l)) Write (6,220)
c     Do loop over grid points 
            Do 20 npt=1,n1
c               Read(8,250,end=10)ref(npt,l,nsp),imf(npt,l,nsp),rhl1,rhl2
c          write(6,400)ref(npt,l,nsp),imf(npt,l,nsp),rhl1,rhl2       
c               If ((ref(npt,l,nsp).eq.rhl2).and.(imf(npt,l,nsp).eq.(-
c     1         rhl1))) GoTo 20
c 10            Write (6,220) nsp,l,npt,ref(npt,l,nsp),rhl1,imf(npt,l,nsp
c     1         ),rhl2
c               STOP
c               last 7 lines were inthe original.Changed mpaez 7may2000
               Read(8,250)ref(npt,l,nsp),imf(npt,l,nsp),rhl1,rhl2
          write(6,400)ref(npt,l,nsp),imf(npt,l,nsp),rhl1,rhl2       
c               If ((ref(npt,l,nsp).eq.rhl2).and.(imf(npt,l,nsp).eq.(-
c     1         rhl1)))
c               GoTo 20
c 10            Write (6,220) nsp,l,npt,ref(npt,l,nsp),rhl1,imf(npt,l,nsp
c     1         ),rhl2
c               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

      Subroutine rcwfn (rho,eta,minl,maxl,fc,fcp,gc,gcp,accur,step)
c *** see RH Landau, Program lpott, 1981
      Implicit  Real*8   (a-h,o-z)
      Real*8 k,k1,k2,k3,k4,m1,m2,m3,m4
      Real*4 rrpace
      Dimension fc(1), fcp(1), gc(1), gcp(1)
c *** Coulomb wavefunctions calculated at r = rho by the
c *** Continued-fraction method of steed   minl,maxl are actual l-values
c *** see barnett feng steed and goldfarb computer physics commun 1974
      pace = step
      acc = accur
      If (pace.lt.100.0) pace = 100.0
      If (acc.lt.1.0e-15.or.acc.gt.1.0e-6) acc = 1.0e-6
      r = rho
      ktr = 1
      lmax = maxl
      lmin1 = minl+1
      xll1 = float(minl*lmin1)
      eta2 = eta*eta
      turn = eta+Sqrt(eta2+xll1)
      If (r.lt.turn.and.abs(eta).ge.1.0e-6) ktr = -1
      ktrp = ktr
      GoTo 20
 10   r = turn
      tf = f
      tfp = fp
      lmax = minl
      ktrp = 1
 20   etar = eta*r
      rho2 = r*r
      pl = float(lmax+1)
      pmx = pl+0.5
c *** Continued fraction for fp(maxl)/f(maxl)  xl is f  xlprime is fp **
      fp = eta/pl+pl/r
      dk = etar*2.0
      del = 0.0
      d = 0.0
      f = 1.0
      k = (pl*pl-pl+etar)*(2.0*pl-1.0)
      If (pl*pl+pl+etar.ne.0.0) GoTo 30
      r = r+1.0e-6
      GoTo 20
 30   h = (pl*pl+eta2)*(1.0-pl*pl)*rho2
      k = k+dk+pl*pl*6.0
      d = 1.0/(d*h+k)
      del = del*(d*k-1.0)
      If (pl.lt.pmx) del = -r*(pl*pl+eta2)*(pl+1.0)*d/pl
      pl = pl+1.0
      fp = fp+del
      If (d.lt.0.0) f = -f
      If (pl.gt.20000.) GoTo 110
      If (abs(del/fp).ge.acc) GoTo 30
      fp = f*fp
      If (lmax.eq.minl) GoTo 50
      fc(lmax+1) = f
      fcp(lmax+1) = fp
c *** Downward recursion to minl for f and fp, arrays gc,gcp are storage
      l = lmax
      Do 40 lp=lmin1,lmax
         pl = float(l)
         gc(l+1) = eta/pl+pl/r
         gcp(l+1) = Sqrt(eta2+pl*pl)/pl
         fc(l) = (gc(l+1)*fc(l+1)+fcp(l+1))/gcp(l+1)
         fcp(l) = gc(l+1)*fc(l)-gcp(l+1)*fc(l+1)
         l = l-1
 40   Continue
      f = fc(lmin1)
      fp = fcp(lmin1)
 50   If (ktrp.eq.-1) GoTo 10
c *** repeat for r = turn If rho lt turn
c *** now obtain p + i.q for minl from Continued fraction (32)
c *** Real arithmetic to facilitate conversion to ibm using r*8
      p = 0.0
      q = r-eta
      pl = 0.0
      ar = -(eta2+xll1)
      ai = eta
      br = 2.0*q
      bi = 2.0
      wi = 2.0*eta
      dr = br/(br*br+bi*bi)
      di = -bi/(br*br+bi*bi)
      dp = -(ar*di+ai*dr)
      dq = (ar*dr-ai*di)
 60   p = p+dp
      q = q+dq
      pl = pl+2.0
      ar = ar+pl
      ai = ai+wi
      bi = bi+2.0
      d = ar*dr-ai*di+br
      di = ai*dr+ar*di+bi
      t = 1.0/(d*d+di*di)
      dr = t*d
      di = -t*di
      h = br*dr-bi*di-1.0
      k = bi*dr+br*di
      t = dp*h-dq*k
      dq = dp*k+dq*h
      dp = t
      If (pl.gt.46000.) GoTo 110
      If (abs(dp)+abs(dq).ge.(abs(p)+abs(q))*acc) GoTo 60
      p = p/r
      q = q/r
c *** solve for fp,g,gp and normalise f  at l=minl
      g = (fp-p*f)/q
      gp = p*g-q*f
      w = 1.0/Sqrt(fp*g-f*gp)
      g = w*g
      gp = w*gp
      If (ktr.eq.1) GoTo 80
      f = tf
      fp = tfp
      lmax = maxl
c *** runge-kutta integration of g(minl) and gp(minl) inwards from turn
c ***             see fox and mayers 1968 pg 202
      If (rho.lt.0.2*turn) pace = 999.0
      r3 = 1.0/3.0d0
      h = (rho-turn)/(pace+1.0)
      h2 = 0.5*h
      rrpace = pace+0.001
      i2 = Ifix(rrpace)
      etah = eta*h
      h2ll = h2*xll1
      s = (etah+h2ll/r)/r-h2
 70   rh2 = r+h2
      t = (etah+h2ll/rh2)/rh2-h2
      k1 = h2*gp
      m1 = s*g
      k2 = h2*(gp+m1)
      m2 = t*(g+m1)
      k3 = h*(gp+m2)
      m3 = t*(g+k2)
      m3 = m3+m3
      k4 = h2*(gp+m3)
      rh = r+h
      s = (etah+h2ll/rh)/rh-h2
      m4 = s*(g+k3)
      g = g+(k1+k2+k2+k3+k4)*r3
      gp = gp+(m1+m2+m2+m3+m4)*r3
      r = rh
      i2 = i2-1
      If (abs(gp).gt.1.0e+36) GoTo 110
c     nb changed exponent to 36 for vax
      If (i2.ge.0) GoTo 70
      w = 1.0/(fp*g-f*gp)
c *** upward recursion from gc(minl) and gcp(minl),stored values are r,s
c *** renormalise fc,fcp for each l-value
 80   gc(lmin1) = g
      gcp(lmin1) = gp
      If (lmax.eq.minl) GoTo 100
      Do 90 l=lmin1,lmax
         t = gc(l+1)
         gc(l+1) = (gc(l)*gc(l+1)-gcp(l))/gcp(l+1)
         gcp(l+1) = gc(l)*gcp(l+1)-gc(l+1)*t
         fc(l+1) = w*fc(l+1)
         fcp(l+1) = w*fcp(l+1)
 90   Continue
      fc(lmin1) = fc(lmin1)*w
      fcp(lmin1) = fcp(lmin1)*w
      Return
 100  fc(lmin1) = w*f
      fcp(lmin1) = w*fp
      Return
 110  w = 0.0
      g = 0.0
      gp = 0.0
      GoTo 80
      End
      Subroutine spbesl (lmaxa,xx,fj)
c *** see RH Landau, Program lpott, 1981
c     spbesl calculates spherical bessel fuctions of 1st      kind for a
c     for all order from zero to lmaxa, at argumrnt xx
      Implicit Real*8  (a-h,o-z)
      Dimension fj(50)
      lmax = lmaxa
      r = xx
      x = xx
      lx = 2.*r
      If (49..lt.lx) lx = 49.
c     calculate small argument limit
      If (r.gt..01) GoTo 10
      fj0 = 1.-x**2/6.+x**4/120.-x**6/42./120.
      If (lmax.gt.1) GoTo 90
      fj(1) = fj0
      fj(2) = x*(1.-x**2/10.+x**4/280.-x**6/30./56./9.)/3.
      GoTo 130
c     regular value of x, start with large l,iterate Down,renormalize
 10   fj0 = sin(x)/x
 20   l2 = max0(lmax+5,lx)
      xl2 = 2.*l2
      fj(l2+1) = 1.e-10
      fj(l2) = 1.e-10*((xl2+1.)/x-x/(xl2+3.))
      l3 = l2-1
      Do 40 ll=1,l3
         l1 = l2-ll
         fl1 = l1
         fj(l1) = (2.*fl1+1.)*fj(l1+1)/x-fj(l1+2)
         rhl = abs(fj(l1))
         If (rhl.le.1.e30) GoTo 40
c     fj too large, scale Down and/or set too small fj*s to 0
         Do 30 l4=l1,l2
            If ((abs(fj(l4)).lt.1.e-26).and.(l4.ne.1)) fj(l4) = 0.
c     nb changed 50 to 26 for vax
            fj(l4) = 1.e-10*fj(l4)
 30      Continue
 40   Continue
c     applay renormaliz factor
      zzj = fj0/fj(1)
      Do 60 l1=1,2
         If (fj(l1).eq.0.) GoTo 50
         If (log10(abs(zzj))+log10(abs(fj(l1))).le.-60.) fj(l1) = 0.
 50      fj(l1) = zzj*fj(l1)
 60   Continue
      If (lmax.le.1) GoTo 130
      l2 = l2-4
      Do 80 l1=3,l2
         If (fj(l1).eq.0.) GoTo 70
         If (log10(abs(zzj))+log10(abs(fj(l1))).le.-60.) fj(l1) = 0.
 70      fj(l1) = zzj*fj(l1)
 80   Continue
      GoTo 130
 90   If (r) 120,100,20
c     x=0 values
 100  Write (6,140)
      l2 = lmax+1
      Do 110 ll=2,l2
         fj(ll) = 0.
 110  Continue
      fj(1) = fj0
      GoTo 130
 120  Write (6,150)
 130  Continue
      Return
c
 140  format (26h spbesl(x) called for x=0,)
 150  format (31h error in spbesl, neg. argument)
      End

      Subroutine legpol (xx,y,n)
      Implicit  Real*8   (a-h,o-z)
c *** see RH Landau, Program lpott, 1981
      Dimension y(50)
      y(1) = 1.
      y(2) = xx
      If (n.gt.2) GoTo 10
      Return
 10   Do 20 i=2,n
         g = xx*y(i)
         y(i+1) = g-y(i-1)+g-(g-y(i-1))/i
 20   Continue
      Return
      End