      SUBROUTINE WAVFN (NSPIN,LMAX,RMEV)
C *** SEE R H 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) GO TO 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))) GO TO 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) GO TO 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) GO TO 120
               IF (NPT.NE.1) GO TO 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
               GO TO 130
C      USE UNDISTORED WABES FOR L GT LSTORE
 120           IF (NPT.NE.N1) GO TO 130
               REUL(L,NSP) = FJ(L)
               AIMUL(L,NSP) = 0.
               IF (ETA.EQ.0.) GO TO 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) GO TO 170
C      L GT LSORE, USE UNDISTORED WAVES
            ZCOS = 1.
            ZSIN = 0.
            ZIDEL = 1.
            GO TO 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
