C
C ERROR CHCKE FOR EIGENSTATE
C CRANCK-NICOLSON SCHEME
C
C (C) TOSHIAKI IITAKA 1994
C REFERENCE: T.Iitaka, Phys. Rev. E49 (1994) 4684.
C
      PROGRAM TSTCN
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NMAX=100)
      PARAMETER(PI=3.14159265358979323)
      COMPLEX*16 P,PH,A1,B1,A2,B2,O,E,F,IU,OVR
      COMPLEX*16 PH0(NMAX)
      DIMENSION P(NMAX),HA(NMAX),HB(NMAX),A1(NMAX),B1(NMAX)
      DIMENSION A2(NMAX),B2(NMAX),O(NMAX),E(NMAX),F(NMAX)
C
C     INITIAL DISPLAY
C
      OPEN(20,FILE='tstcn.dat')
      WRITE(*,*) 'M, ALPHA'
      READ(*,*) M,ALPHA
      IU=(0.0,1.0)
      DX=1
      DX2=DX*DX
      EMAX=1/DX2
      DT=ALPHA/(EMAX*COS(M*PI/(NMAX+1)))
      MTIME=(2*PI/ALPHA) * 10
      MTIME=MIN(MTIME,10000)
C
C     CALCULATION OF HAMILTONIAN
C
      DO 50 N=1,NMAX
C        HA(N) = V(X)+1/DX2 - E0
C        HA(N) = 1/DX2-E0
        HA(N)=0
        HB(N) = -0.5/DX2
        A1(N) = 1.-0.5*IU*DT*HA(N)
        A2(N) = 1.+0.5*IU*DT*HA(N)
        B1(N) =   -0.5*IU*DT*HB(N)
        B2(N) =   +0.5*IU*DT*HB(N)
   50 CONTINUE
C
C     SETTING INITIAL WAVE FUNCTION
C
      CALL EIGEN(M,DX,PH0,EM)
      DO 60 N=1,NMAX
      P(N)=PH0(N)
   60 CONTINUE
C
C     TIME EVOLUTION LOOP 1000
C
      DO 1000 ITIME=0,MTIME
      T=ITIME*DT
C
C      OUTPUT WAVE FUNCTION
C
      IF(MOD(ITIME,10).EQ. 1) THEN
      ANORM=0
      OVR=0.0
      DO 150 N=1,NMAX
        ANORM=ANORM+ABS(P(N))**2
        OVR=OVR+DCONJG(P(N))*PH0(N)*EXP(-IU*EM*T)
  150   CONTINUE
      WRITE(*,100) ITIME,EMAX*T,ANORM-1, 
     &1.- ABS(OVR),DATAN2(DIMAG(OVR),DBLE(OVR))
     &,ABS((EM*DT)**3/12.) * ITIME
      WRITE(20,100) ITIME,EMAX*T,ABS(ANORM-1),
     &ABS(1- ABS(OVR)),ABS(DATAN2(DIMAG(OVR),DBLE(OVR)))
     &,ABS((EM*DT)**3/12.) * ITIME
  100 FORMAT(1H ,I4,30E15.7)
      ENDIF
C
C     MULTIPLICATION OF SYMETRIC TRIDIAGONAL MATRIX
C                      O=(AB)*P
C
      O(1) = A1(1)*P(1)+B1(1)*P(2)
      DO 70 N=2,NMAX-1
         O(N) = B1(N-1)*P(N-1)+A1(N)*P(N)+B1(N)*P(N+1)
   70 CONTINUE
         O(NMAX) = B1(NMAX-1)*P(NMAX-1)+A1(NMAX)*P(NMAX)
C
C     INVERSION OF SYNMETRIC TRIDIAGONAL MATRIX 
C                      P=O/(AB)
C
      E(1) = -A2(1)/B2(1)
      F(1) =  O(1)/B2(1)
      DO 80 N=2,NMAX
        E(N)= -(B2(N-1)/E(N-1)+A2(N)) / B2(N)
        F(N)= O(N)/B2(N) + (B2(N-1)/B2(N))*(F(N-1)/E(N-1))
   80 CONTINUE
      P(NMAX)=-F(NMAX)/E(NMAX)
      DO 90 N=NMAX-1,1,-1
        P(N)=(P(N+1)-F(N))/E(N)
   90 CONTINUE
 1000 CONTINUE
      STOP
      END
C
C     SETTING INITIAL WAVE FUNCTION
C
      SUBROUTINE EIGEN(M,DX,PH0,EM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(PI=3.14159265358979323)
      PARAMETER(NMAX=100,IMAX=1000)
      PARAMETER(EPS=1E-14)
      COMPLEX*16 PH0(NMAX),IU
      DIMENSION HA(NMAX),HB(NMAX-1),V(NMAX)
      IU=(0.0,1.0)
      DO 60 N=1,NMAX
        PH0(N)= SIN((PI*M*N)/(NMAX+1))
   60 CONTINUE
C     NORMALIZATION OF WAVE FUNCTION
      PA = 0.
      DO 100 N=1,NMAX
        PA = PA + ABS(PH0(N))**2
  100 CONTINUE
      PA = SQRT(PA)
      DO 110 N=1,NMAX
        PH0(N)=PH0(N)/PA
  110 CONTINUE
C      EM=-4*(-0.5/DX/DX) * SIN((PI*M)/(NMAX+1)/2.0)**2
      EM=2*(-0.5/DX/DX) * COS((PI*M)/(NMAX+1))
      RETURN
      END