PROGRAM KAP15

*---------------------------------------------------------------------
*     Program KAP15 solves the radial SCHR™DINGER equation with
*     angular momentum l as a matrix equation of dimension D. As
*     representation basis the harmonic oscillator states are used.
*     The matrix elements of the HAMILTON operator H are first
*     calculated, then the eigenvalue equation is solved and the wave
*     functions are calculated from the eigenvectors. The à-à system
*     serves as example.
*    --Fixed quantities--
*     H2M:    hý/(8ãým) in units of 1MeVfmý
*             (h Planck's constant, m (reduced) mass,
*             in the à-à example H2M=10.375D0)
*             {DOUBLE PRECISION constant}
*     NMAX:   Maximum dimension of the function basis (NMAX=20)
*             {INTEGER constant}
*     B:      Upper integration limit b in units of 1 fm (B=10.D0)
*             {DOUBLE PRECISION constant}
*     IB:     Number of integration steps in the integration
*             interval [0,b] (IB=500)
*             {INTEGER constant}
*    --Input quantities--
*     L:      Angular momentum quantum number l
*             {INTEGER variable}
*     NDIM:   Dimension of the function basis (1 ó NDIM ó NMAX)
*             {INTEGER variable}
*     A:      Oscillator width parameter in units of 1/fmý
*             {DOUBLE PRECISION variable}
*    --Output quantities--
*     HMAT:   Matrix of the HAMILTON operator
*             {DOUBLE PRECISION array of dimensions (NMAX,NMAX)}
*     E:      Energy eigenvalue in units of 1 MeV
*             {DOUBLE PRECISION array of dimension (NMAX)}
*     U:      Radial wave functions u,
*             U(I,N) is the function value of the wave function at
*             the point R(I) pertaining to the eigenvalue E(N)
*             {DOUBLE PRECISION array of dimensions (0:IB,NMAX)}
*    --Other quantities--
*     R:      Distance in units of 1 fm,
*             R(I)=I*H, H=B/IB, I=0,1,2,...,IB
*             {DOUBLE PRECISION array of dimension (0:IMAX)}
*     VR:     Potential V in units of 1 MeV,
*             VR(I) is the potential value at distance R(I)
*             {DOUBLE PRECISION array of dimension (0:IMAX)}
*     VNL:    Radial oscillator function with the radial quantum
*             number n and the angular momentum quantum number l,
*             VNL(I,N) is the function value at the point R(I)
*             {DOUBLE PRECISION array of dimensions (0:IB,NMAX)}
*    --Subprograms called--
*     EIGEN, RADOS, GOPEN, GMOVE, GDRAW, GCLOSE, CURSOR  {SUBROUTINES}
*     V  {DOUBLE PRECISION FUNCTION}
*---------------------------------------------------------------------

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (IB=500, NMAX=20, B=10.D0, H2M=10.375D0)
      DIMENSION E(NMAX), H2(NMAX), H3(NMAX), OS(NMAX)
      DIMENSION HMAT(NMAX,NMAX), UMAT(NMAX,NMAX), H1(NMAX,NMAX)
      DIMENSION U(0:IB,NMAX), VNL(0:IB,NMAX), R(0:IB), VR(0:IB)

*    >>>>>Input/output
      CHARACTER CONT*1
      PARAMETER (ILINE=25, IOMAX=105)
      CALL CURSOR(1)
      WRITE(*,'(T29,A)')         'ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿'
      WRITE(*,'(T29,A)')         '³ C H A P T E R   1 5 ³'
      WRITE(*,'(T29,A/)')        'ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ'
      WRITE(*,'(T3,A/8(T2/))')
     &'SOLUTION OF THE SCHR™DINGER EQUATION IN HARMONIC OSCILLATOR REPRE
     &SENTATION'
      WRITE(*,'(T2,A)')
     &'Please input values for the following quantities:'
      IFLAG=1
1000  CONTINUE
      WRITE(*,'(T2/T2,A)')
     &'Angular momentum quantum number L  (0 ó L ó 10):  L = '
      READ(*,*,IOSTAT=IOS) L
      IF((L.LT.0).OR.(L.GT.10).OR.(IOS.NE.0)) THEN
       WRITE(*,'(T2,A/)') 'Input error !'
       GOTO 1000
      ENDIF
2000  CONTINUE
      WRITE(*,'(T2,A)')       'Oscillator width parameter A in units of
     &1/fmý  (0 < A ó 2.0):  A = '
      READ(*,*,IOSTAT=IOS) A
      IF((A.LE.0).OR.(A.GT.2.D0).OR.(IOS.NE.0)) THEN
       WRITE(*,'(T2,A/)') 'Input error !'
       GOTO 2000
      ENDIF
3000  CONTINUE
      WRITE(*,'(T2,A,I2,A)') 'Dimension NDIM of the representation basis
     &  (1 ó NDIM ó ',NMAX,'):  NDIM = '
      READ(*,*,IOSTAT=IOS) NDIM
      IF((NDIM.LT.1).OR.(NDIM.GT.NMAX).OR.(IOS.NE.0)) THEN
       WRITE(*,'(T2,A/)') 'Input error !'
       GOTO 3000
      ENDIF
*    <<<<<

      H=B/IB
      DO 20 I=0,IB
       R(I)=I*H
       VR(I)=V(R(I))
       CALL RADOS(R(I),A,L,NDIM,OS)
       DO 10 N=1,NDIM
        VNL(I,N)=OS(N)
10     CONTINUE
20    CONTINUE
      HQOM=4.D0*A*H2M
      OM2=8.D0*A*A*H2M
      DO 60 N=1,NDIM
       DO 50 NP=1,N
        IF(NP.EQ.N) THEN
         HMAT(N,NP)=(2.D0*N+L-0.5D0)*HQOM/2.D0
        ELSE IF(NP.EQ.N-1) THEN
          HMAT(N,NP)=SQRT(NP*(NP+L+0.5D0))*HQOM/2.D0
        ELSE
          HMAT(N,NP)=0.D0
        ENDIF
        DO 30 I=1,IB-1,2
         F=VNL(I,N)*VR(I)*VNL(I,NP)
         HMAT(N,NP)=HMAT(N,NP)+4.D0*H/3.D0*F
30      CONTINUE
        DO 40 I=2,IB-1,2
         F=VNL(I,N)*VR(I)*VNL(I,NP)
         HMAT(N,NP)=HMAT(N,NP)+2.D0*H/3.D0*F
40      CONTINUE
        HMAT(NP,N)=HMAT(N,NP)
50     CONTINUE
60    CONTINUE

*    >>>>>Input/output
4000  CONTINUE
      WRITE(*,'(T2/T2,A)')'Output of the matrix HMAT required  (Y/N) ? '
      READ(*,'(A1)') CONT
      IF(INDEX('Yy ',CONT).NE.0) THEN
       CALL CURSOR(1)
       WRITE(*,'(T2,A,I2,A,F6.3,A,I2,A/)') 'Matrix elements of HMAT (rec
     &all: L=',L,', A=',A,', NDIM=',NDIM,'):'
       NLOOP=0
5000   CONTINUE
       NLOOP=NLOOP+10
       NNLOOP=NLOOP-9
       IF(NLOOP.GT.N) NLOOP=NDIM
       NR=0
6000   CONTINUE
       NR=NR+5
       NNR=NR-4
       IF(NR.GT.NDIM) NR=NDIM
       WRITE(*,'(T16,(5(A,I2,A)))')
     & (' Column ',NP,'   ',NP=NNR,NR)
       WRITE(*,'(T2)')
       DO 7000 N=NNLOOP,NLOOP
        WRITE(*,'(T2,A,I2,A,5F13.6)')
     &  '  Row ',N,'   ',(HMAT(N,NP),NP=NNR,NR)
7000   CONTINUE
       IF(NR.LT.NDIM) THEN
        WRITE(*,'(T2/T2,A)') 'Continue ?  Please press ENTER key ...'
        READ(*,'(A1)') CONT
        CALL CURSOR(3)
        GOTO 6000
       ELSEIF(NLOOP.LT.NDIM) THEN
        WRITE(*,'(T2/T2,A)') 'Continue ?  Please press ENTER key ...'
        READ(*,'(A1)') CONT
        CALL CURSOR(3)
        GOTO 5000
       ENDIF
8000   CONTINUE
       WRITE(*,'(T2/T2,A,I2,A/T2,A)')
     & 'Rerun calculation of the matrix HMAT for L = ',L,' with',
     & 'different values for A and NDIM  (Y/N) ? '
       READ(*,'(A1)') CONT
       IF(INDEX('Yy ',CONT).NE.0) THEN
        GOTO 2000
       ELSEIF(INDEX('Nn',CONT).EQ.0) THEN
        GOTO 8000
       ENDIF
      ELSEIF(INDEX('Nn',CONT).EQ.0) THEN
       GOTO 4000
      ENDIF
*    <<<<<

      IF(NDIM.GT.3) WRITE(*,'(T2,A)') 'Please wait ... '
      CALL EIGEN(HMAT,NMAX,NDIM,E,UMAT,H1,H2,H3)
      DO 90 NP=1,NDIM
       DO 80 I=0,IB
        U(I,NP)=0.D0
        DO 70 N=1,NDIM
         U(I,NP)=U(I,NP)+UMAT(N,NP)*VNL(I,N)
70      CONTINUE
80     CONTINUE
90    CONTINUE

*    >>>>>Input/output
      CALL CURSOR(1)
      WRITE(*,'(T2,A,I2,A,F6.3,A,I2,A/)')        'Energy eigenvalues (re
     &call: L=',L,', A=',A,', NDIM=',NDIM,'):'
      WRITE(*,'(T2,A/)')
     &'No.       Eigenvalue                   No.       Eigenvalue'
      DO 9000 NR=1,NDIM/2
       NNR=NR+NDIM-(NDIM/2)
       WRITE(*,'(T2,I2,F18.6,A,T41,I2,F18.6,A)')
     & NR,E(NR),' MeV',NNR,E(NNR),' MeV'
9000  CONTINUE
      IF(NDIM.NE.2*(NDIM/2))
     & WRITE(*,'(T2,I2,F18.6,A)') NDIM-(NDIM/2),E(NDIM-(NDIM/2)),' MeV'
10000 CONTINUE
      WRITE(*,'(T2/T2,A)')
     &'Output of radial wave functions required  (Y/N) ? '
      READ(*,'(A1)') CONT
      IF(INDEX('Yy ',CONT).NE.0) THEN
       RMAX=MIN(B,SQRT((2.D0*NDIM+1.5D0+L)/A)+3.D0)
11000  CONTINUE
       WRITE(*,'(T2,A)')
     & 'For which eigenvalue ? Input number of eigenvalue --> '
       READ(*,*,IOSTAT=IOS) NR
       IF((NR.LT.1).OR.(NR.GT.NDIM).OR.(IOS.NE.0)) THEN
        WRITE(*,'(T2,A/)') 'Input error !'
        GOTO 11000
       ENDIF
       UMIN=U(0,NR)
       UMAX=U(0,NR)
       DO 12000 I=1,IB
        UMIN=MIN(UMIN,U(I,NR))
        UMAX=MAX(UMAX,U(I,NR))
12000  CONTINUE
       CALL GOPEN(0.D0,RMAX,-1.D0,UMIN,UMAX,-1.D0,
     &  'Distance [fm]','Radial wave function',0,1)
       IF(IFLAG.EQ.1) THEN
        CALL GCOLOR('Please choose colour for the wave function',ICOLOR)
        IFLAG=0
       ENDIF
       CALL GMOVE(1,R(0),U(0,NR),ICOLOR,1)
       DO 13000 I=1,IB
        CALL GDRAW(1,R(I),U(I,NR))
13000  CONTINUE
       CALL GCLOSE
14000  CONTINUE
       WRITE(*,'(T2,A)') 'Numerical output required  (Y/N) ? '
       READ(*,'(A1)') CONT
       IF(INDEX('Yy ',CONT).NE.0) THEN
        IS=(IB-1)/IOMAX+1
        JLINE=0
        DO 15000 I=0,IB,IS
         IF(JLINE.EQ.0) WRITE(*,'(T14,A,T33,A/)')
     &    'Distance [fm]','Radial wave function'
         WRITE(*,'(T2,F20.2,F25.6)') R(I),U(I,NR)
         JLINE=JLINE+1
         IF((JLINE.EQ.ILINE-4).OR.(I.GT.IB-IS)) THEN
          WRITE(*,'(T2/T2,A)') 'Continue ?  Please press ENTER key ...'
          READ(*,'(A1)') CONT
          JLINE=0
          CALL CURSOR(1)
         ENDIF
15000   CONTINUE
       ELSEIF(INDEX('Nn',CONT).EQ.0) THEN
        GOTO 14000
       ENDIF
16000  CONTINUE
       WRITE(*,'(T2,A)')
     & 'Output of further wave functions required  (Y/N) ? '
       READ(*,'(A1)') CONT
       IF(INDEX('Yy ',CONT).NE.0) THEN
        GOTO 11000
       ELSEIF(INDEX('Nn',CONT).EQ.0) THEN
        GOTO 16000
       ENDIF
      ELSEIF(INDEX('Nn',CONT).EQ.0) THEN
       GOTO 10000
      ENDIF
17000 CONTINUE
      WRITE(*,'(T2/T2,A,I2,A/T2,A)')
     &'Rerun program for L = ',L,' with',
     &'different values for A and NDIM  (Y/N) ? '
      READ(*,'(A1)') CONT
      IF(INDEX('Yy ',CONT).NE.0) THEN
       GOTO 2000
      ELSEIF(INDEX('Nn',CONT).EQ.0) THEN
       GOTO 17000
      ENDIF
18000 CONTINUE
      WRITE(*,'(T2/T2,A)')
     &'Rerun program with completely new values  (Y/N) ? '
      READ(*,'(A1)') CONT
      IF(INDEX('Yy ',CONT).NE.0) THEN
       GOTO 1000
      ELSEIF(INDEX('Nn',CONT).EQ.0) THEN
       GOTO 18000
      ENDIF
      WRITE(*,'(T2//T2,A///)') 'Program ends !'
*    <<<<<

      END


      DOUBLE PRECISION FUNCTION V(R)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (V0=122.694D0, BETA=0.22D0)
      PARAMETER (EXPMAX=50.D0)
      IF(BETA*R*R.LE.EXPMAX) THEN
       V=-V0*EXP(-BETA*R*R)
      ELSE
       V=0.D0
      ENDIF
      END


      SUBROUTINE EIGEN(A,NMAX,N,ALPHA,X,H1,H2,H3)
*---------------------------------------------------------------------
*     The subroutine EIGEN calculates the eigenvalues and eigenvectors
*     of a real, symmetric matrix A for the eigenvalue problem Ax=àx
*     using the JACOBI method.
*    --Input parameters--
*     A:      Real, symmetric matrix
*             {DOUBLE PRECISION array of dimensions (NMAX,NMAX)}
*     NMAX:   Maximum dimension
*             {INTEGER constant}
*     N:      Dimension of the matrix A (1 ó N ó NMAX)
*             {INTEGER variable}
*     H1:     Auxiliary array
*             {DOUBLE PRECISION array of dimensions (NMAX,NMAX)}
*     H2:     Auxiliary array
*             {DOUBLE PRECISION array of dimension (NMAX)}
*     H3:     Auxiliary array
*             {DOUBLE PRECISION array of dimension (NMAX)}
*    --Output parameters-
*     ALPHA:  Eigenvalues à, where ALPHA(1) ó ... ó ALPHA(N)
*             {DOUBLE PRECISION array of dimension (NMAX)}
*     X:      Normalised eigenvectors x, where (X(1,K),..,X(N,K))
*             is the eigenvector for the eigenvalue ALPHA(K)
*             {DOUBLE PRECISION array of dimensions (NMAX,NMAX)}
*---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (LIMIT=50, EPS=1.D-14)
      DIMENSION ALPHA(NMAX), H2(NMAX), H3(NMAX)
      DIMENSION A(NMAX,NMAX), X(NMAX,NMAX), H1(NMAX,NMAX)
      IF((N.LT.1).OR.(N.GT.NMAX)) STOP
     &'Error in EIGEN: 1 ó N ó NMAX not satisfied'
      DO 20 N1=1,N
       DO 10 N2=1,N
        H1(N1,N2)=A(N1,N2)
        IF(N1.EQ.N2) THEN
         X(N1,N2)=1.D0
        ELSE
         X(N1,N2)=0.D0
        ENDIF
10     CONTINUE
20    CONTINUE
      DO 30 N1=1,N
       ALPHA(N1)=H1(N1,N1)
       H2(N1)=ALPHA(N1)
       H3(N1)=0.D0
30    CONTINUE
      IF(N.GT.1) THEN
       DO 160 I=1,LIMIT
        SM=0.D0
        DO 50 N1=1,N-1
         DO 40 N2=N1+1,N
          SM=SM+ABS(H1(N1,N2))
40       CONTINUE
50      CONTINUE
        IF(SM.GT.EPS) THEN
         IF(I.LT.4) THEN
          TRESH=SM/(5.D0*N*N)
         ELSE
          TRESH=0.D0
         ENDIF
         DO 110 N1=1,N-1
          DO 100 N2=N1+1,N
           G=100.D0*ABS(H1(N1,N2))
           IF((I.GT.4).AND.(ABS(ALPHA(N1))+G.LT.ABS(ALPHA(N1))+EPS)
     &      .AND.(ABS(ALPHA(N2))+G.LT.ABS(ALPHA(N2))+EPS)) THEN
            H1(N1,N2)=0.D0
           ELSEIF(ABS(H1(N1,N2)).GT.TRESH) THEN
            H=ALPHA(N2)-ALPHA(N1)
            IF(ABS(H)+G.LT.ABS(H)+EPS) THEN
             T=H1(N1,N2)/H
            ELSE
             THETA=H/(2.D0*H1(N1,N2))
             T=1.D0/(ABS(THETA)+SQRT(1.D0+THETA*THETA))
             IF(THETA.LT.0.D0) T=-T
            ENDIF
            C=1.D0/SQRT(1.D0+T*T)
            S=T*C
            TAU=S/(1.D0+C)
            H=T*H1(N1,N2)
            H3(N1)=H3(N1)-H
            H3(N2)=H3(N2)+H
            ALPHA(N1)=ALPHA(N1)-H
            ALPHA(N2)=ALPHA(N2)+H
            H1(N1,N2)=0.D0
            IF(N1.GT.1) THEN
             DO 60 J=1,N1-1
              G=H1(J,N1)
              H=H1(J,N2)
              H1(J,N1)=G-S*(H+G*TAU)
              H1(J,N2)=H+S*(G-H*TAU)
60           CONTINUE
            ENDIF
            IF((N1+1).LE.(N2-1)) THEN
             DO 70 J=N1+1,N2-1
              G=H1(N1,J)
              H=H1(J,N2)
              H1(N1,J)=G-S*(H+G*TAU)
              H1(J,N2)=H+S*(G-H*TAU)
70           CONTINUE
            ENDIF
            IF(N2.LT.N) THEN
             DO 80 J=N2+1,N
              G=H1(N1,J)
              H=H1(N2,J)
              H1(N1,J)=G-S*(H+G*TAU)
              H1(N2,J)=H+S*(G-H*TAU)
80           CONTINUE
            ENDIF
            DO 90 J=1,N
             G=X(J,N1)
             H=X(J,N2)
             X(J,N1)=G-S*(H+G*TAU)
             X(J,N2)=H+S*(G-H*TAU)
90          CONTINUE
           ENDIF
100       CONTINUE
110      CONTINUE
         DO 120 N1=1,N
          H2(N1)=H2(N1)+H3(N1)
          ALPHA(N1)=H2(N1)
          H3(N1)=0.D0
120      CONTINUE
        ELSE
         DO 150 K=1,N-1
          DO 140 L=K+1,N
           IF(ALPHA(K).GT.ALPHA(L)) THEN
            AUX=ALPHA(K)
            ALPHA(K)=ALPHA(L)
            ALPHA(L)=AUX
            DO 130 J=1,N
             H3(J)=X(J,K)
             X(J,K)=X(J,L)
             X(J,L)=H3(J)
130         CONTINUE
           ENDIF
140       CONTINUE
150      CONTINUE
         GOTO 170
        ENDIF
160    CONTINUE
       STOP 'Error in EIGEN: no convergence'
      ENDIF
170   CONTINUE
      END


      SUBROUTINE RADOS(R,A,L,NMAX,OS)
*---------------------------------------------------------------------
*     Subroutine RADOS calculates the radial oscillator functions for
*     the angular momentum quantum number l and all radial quantum
*     numbers from 1 to nmax at the point r.
*    --Input parameters--
*     R:      Distance r (R ò 0.D0)
*             {DOUBLE PRECISION variable}
*     A:      Oscillator width parameter
*             {DOUBLE PRECISION variable}
*     L:      Angular momentum quantum number l (L ò 0)
*             {INTEGER variable}
*     NMAX:   Maximum radial quantum number nmax (NMAX ò 1)
*             {INTEGER constant}
*    --Output parameter--
*     OS:     Radial oscillator function,
*             OS(N) is the function value at the point r for the
*             radial quantum number n
*             {DOUBLE PRECISION array of dimension (NMAX)}
*---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (PI=3.1415926535D0)
      PARAMETER (EXPMAX=700.D0)
      DIMENSION OS(NMAX)
      IF((R.LT.0.D0).OR.(L.LT.0).OR.(NMAX.LT.1)) STOP
     &'Error in RADOS: False value for R, L or NMAX !'
      R0=1.D0/SQRT(2.D0*A)
      X=(R*R)/(R0*R0)
      CL=4.D0/(SQRT(PI)*R0)*X
      IF(L.GE.1) THEN
       DO 10 J=1,L
        CL=2*CL/(2*J+1)*X
10     CONTINUE
      ENDIF
      IF(X/2.LE.EXPMAX) THEN
       OS(1)=SQRT(CL)*EXP(-X/2)
      ELSE
       OS(1)=0.D0
      ENDIF
      IF(NMAX.GE.2) OS(2)=SQRT(2.D0/(2*L+3))*OS(1)*(L+1.5D0-X)
      IF(NMAX.GE.3) THEN
       DO 20 N=2,NMAX-1
        OS(N+1)=SQRT(1.D0/(N*(N+L+0.5D0)))*
     &   ( (2*N+L-0.5D0-X)*OS(N) - SQRT((N-1)*(N+L-0.5D0))*OS(N-1) )
20     CONTINUE
      ENDIF
      END