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 SCCN
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NMAX=100)
      PARAMETER(PI=3.14159265358979323)
      COMPLEX*16 P,A1,B1,A2,B2,O,E,F,IU,OVR
      COMPLEX*16 PH0(NMAX),PH(NMAX)
      DIMENSION P(NMAX),HA(NMAX),HB(NMAX),A1(NMAX),B1(NMAX),V(NMAX)
      DIMENSION A2(NMAX),B2(NMAX),O(NMAX),E(NMAX),F(NMAX)
      DIMENSION HAN(NMAX),HBN(NMAX)
C
C     INITIAL DISPLAY
C
      OPEN(20,FILE='sccn.dat')
      WRITE(*,*) 'ALPHA'
      READ(*,*) ALPHA
      IU=(0.0,1.0)
      DX=1
      DX2=DX*DX
      EMAX=1/DX2+1
      DT=ALPHA/EMAX
      MTIME=(2*PI/ALPHA) * 40
      MTIME=MIN(MTIME,100000)
C
C     CALCULATION OF HAMILTONIAN
C
      DO 50 N=1,NMAX
      IF(45.LE.N .AND. N.LE.55) THEN
        V(N)=0.05*COS((I-50)/5.0 * PI/2)**2
c        V(N)=0.05
      ELSE
        V(N)=0
      ENDIF
       HA(N) = V(N) 
c       HA(N) = V(N) + 1/DX2
C      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 NORMALIZATION OF HAMILTONIAN
C
      VMAX=0
      VMIN=0
      DO 10 N=1,NMAX
      IF(V(N).GT.VMAX) VMAX=V(N)
      IF(V(N).LT.VMIN) VMIN=V(N)
   10 CONTINUE
C      AKMAX= 0.5 * (PI/DX)**2
C      EMIN=VMIN
C      EMAX=VMAX+AKMAX
      VMAX=1.0
      WRITE(*,*) 'VMAX=',VMAX
      EMAX=1/DX2 +VMAX
      EMIN=-1/DX2
      EGRID=EMAX-EMIN
      DO 20 N=1,NMAX
      HAN(N)=(HA(N)-EGRID/2-EMIN)*2/EGRID
      IF(N.LT.NMAX) HBN(N)=HB(N)*2/EGRID
   20 CONTINUE
      WRITE(*,*) 'EGRID=',EGRID
      WRITE(*,*) 'TIME,M'
C
C SET INITAIL FUNCTION
C
C
C     SETTING INITIAL WAVE FUNCTION
C
C      WRITE(*,*)'INPUT THE CENTER OF WAVE PACKET 
C     &           AND ITS WIDTH AND MOMENTUM'
C      READ(*,*) X0,SG,PA
      X0=0.25*DX*NMAX
      SG=0.1*DX*NMAX
      PA=0.1*2*PI/DX
      DO 60 N=1,NMAX
        X=N*DX
        PH0(N)= EXP(-0.5*(X-X0)**2/SG**2) * EXP(IU*PA*(X-X0))
   60 CONTINUE
C
C     NORMALIZATION OF WAVE FUNCTION
C
      PA = 0.
      DO 101 N=1,NMAX
        PA = PA + ABS(PH0(N))**2
  101 CONTINUE
      PA = SQRT(PA)
      DO 110 N=1,NMAX
        PH0(N)=PH0(N)/PA
  110 CONTINUE
      DO 81 N=1,NMAX
        P(N)=PH0(N)
   81 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,30).EQ. 0) THEN
      CALL CH(HAN,HBN,EGRID,EMIN,PH0,T,PH)
      ANORM=0
      OVR=0.0
      DO 150 N=1,NMAX
        ANORM=ANORM+ABS(P(N))**2
        OVR=OVR+DCONJG(P(N))*PH(N)
  150   CONTINUE
      WRITE(*,100) ITIME,EMAX*T,ANORM-1, 
     &1.- ABS(OVR),DATAN2(DIMAG(OVR),DBLE(OVR))
      WRITE(20,100) ITIME,EMAX*T,ABS(ANORM-1),
     &ABS(1- ABS(OVR)),ABS(DATAN2(DIMAG(OVR),DBLE(OVR)))
  100 FORMAT(1H ,I7,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
C      WRITE(*,*) B2(N),E(N-1)
        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 SUBROUTINE CH
C
      SUBROUTINE CH(HA,HB,EGRID,EMIN,PH0,T,PH)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(PI=3.14159265358979323)
      PARAMETER(NMAX=100,IMAX=1000)
      PARAMETER(EPS=1E-16)
      COMPLEX*16 P(NMAX),Q(NMAX),R(NMAX),PH(NMAX),PH0(NMAX),IU
      DIMENSION HA(NMAX),HB(NMAX-1),V(NMAX)
      IU=(0.0,1.0)
C
      DO 20 N=1,NMAX
      P(N)=PH0(N)
   20 CONTINUE
      Q(1)   = -IU*(HA(1)*P(1)+HB(1)*P(2))
      DO 30 N=2,NMAX-1
      Q(N)   = -IU*(HB(N-1)*P(N-1)+HA(N)*P(N)+HB(N)*P(N+1))
   30 CONTINUE
      Q(NMAX)= -IU*(HB(NMAX-1)*P(NMAX-1)+HA(NMAX)*P(NMAX))
C
      AN=DJBES(0,EGRID*T/2)
      DO 40 N=1,NMAX
      PH(N)=AN*P(N)
   40 CONTINUE
      AN=2*DJBES(1,EGRID*T/2)
      DO 45 N=1,NMAX
      PH(N)=PH(N)+AN*Q(N)
   45 CONTINUE
      DO 1000 I=2, IMAX, 3
C
      AN=2*DJBES(I,EGRID*T/2)
      R(1) = -2*IU*(HA(1)*Q(1)+HB(1)*Q(2))+P(1)
      PH(1)=PH(1) +AN*R(1)
      DO 70 N=2,NMAX-1
      R(N)=-2*IU*(HB(N-1)*Q(N-1)+HA(N)*Q(N)+HB(N)*Q(N+1))+P(N)
      PH(N)=PH(N)+AN*R(N)
   70 CONTINUE
      R(NMAX)=-2*IU*(HB(NMAX-1)*Q(NMAX-1)+HA(NMAX)*Q(NMAX))+P(NMAX)
      PH(NMAX)=PH(NMAX)+AN*R(NMAX)
C
      AN=2*DJBES(I+1,EGRID*T/2)
      P(1) = -2*IU*(HA(1)*R(1)+HB(1)*R(2))+Q(1)
      PH(1)=PH(1) +AN*P(1)
      DO 71 N=2,NMAX-1
      P(N)=-2*IU*(HB(N-1)*R(N-1)+HA(N)*R(N)+HB(N)*R(N+1))+Q(N)
      PH(N)=PH(N)+AN*P(N)
   71 CONTINUE
      P(NMAX)=-2*IU*(HB(NMAX-1)*R(NMAX-1)+HA(NMAX)*R(NMAX))+Q(NMAX)
      PH(NMAX)=PH(NMAX)+AN*P(NMAX)
C
      AN=2*DJBES(I+2,EGRID*T/2)
      Q(1) = -2*IU*(HA(1)*P(1)+HB(1)*P(2))+R(1)
      PH(1)=PH(1) +AN*Q(1)
      DO 72 N=2,NMAX-1
      Q(N)=-2*IU*(HB(N-1)*P(N-1)+HA(N)*P(N)+HB(N)*P(N+1))+R(N)
      PH(N)=PH(N)+AN*Q(N)
   72 CONTINUE
      Q(NMAX)=-2*IU*(HB(NMAX-1)*P(NMAX-1)+HA(NMAX)*P(NMAX))+R(NMAX)
      PH(NMAX)=PH(NMAX)+AN*Q(NMAX)
C      WRITE(*,*) 'I=',I, AN
      IF(ABS(AN).LT.EPS) GOTO 2000
 1000 CONTINUE
      WRITE(*,*) 'CHEBYSHEV NOT CONVERGE'
 2000 CONTINUE
      DO 80 N=1,NMAX
      PH(N)=PH(N)*EXP(-IU*(EGRID/2+EMIN)*T)
   80 CONTINUE
      RETURN
      END
C
C  BESSEL FUNCTION DJBES(N,X)
C
      REAL FUNCTION DJBES * 8 (N,X)                                     
      REAL *8 X,Z,BJ,QJ,S,T1,T2,T3,ONE,D55                              
C      CALL OVERFL(KJ)                                                   
      ONE=1.0D0                                                         
      D55=1.0D-55                                                       
      NN=IABS(N)                                                        
      XA=DABS(X)                                                        
      IF(NN-30000)10,900,900                                            
   10 IF(XA-10.0)11,13,13                                               
   11 IF(XA-1.0)12,12,15                                                
   12 IF(XA-0.00002)1,1,16                                              
   13 IF(XA-100.0)17,14,14                                              
   14 IF(XA-30000.0)18,900,900                                          
    1 T1=0.5D0*X                                                        
      IF(NN) 3,2,3                                                      
    2 DJBES = ONE-T1*T1                                                 
      GO TO 1000                                                        
    3 IF(XA .LE. 1.0D-77) GO TO 500                                     
      T2=ONE                                                            
      T3=ONE                                                            
      DO 5 I=1,NN                                                       
      IF(DABS(T3) .LE. 1.0D-77*DABS(T2/T1)) GO TO 500                   
      T3=T3*T1/T2                                                       
      T2=T2+ONE                                                         
    5 CONTINUE                                                          
      BJ=T3*(ONE-T1*T1/T2)                                              
      GO TO 300                                                         
   15 L=1.4*XA+14.0                                                     
      GO TO 20                                                          
   16 L=14                                                              
      GO TO 20                                                          
   17 L=0.27*XA+27.0                                                    
      GO TO 20                                                          
   18 L=0.073*XA+47.0                                                   
   20 Z=2.0D0/X                                                         
      NM=MAX0(NN,IFIX(XA))+L                                            
      T3=0.0D0                                                          
      T2=1.0D-75                                                        
      S=0.0D0                                                           
      IF( MOD(NM,2)) 22,21,22                                           
   21 NM=NM+1                                                           
   22 DO 100 I=1,NM,2                                                   
      K=NM-I+1                                                          
      T1=DFLOAT(K+1)*T2*Z-T3                                            
      IF(NN-K) 40,30,40                                                 
   30 QJ=T1                                                             
   40 K=K-1                                                             
      T3=T2                                                             
      T2=T1                                                             
      T1=DFLOAT(K+1)*T2*Z-T3                                            
      IF(NN-K) 50,45,50                                                 
   45 QJ=T1                                                             
   50 S=S+T1                                                            
      IF( DABS(S)-1.0D55 ) 80,60,60                                     
   60 T1=T1*D55                                                         
      T2=T2*D55                                                         
      S =S*D55                                                          
      IF (NN-K) 80,70,70                                                
   70 QJ=QJ*D55                                                         
   80 T3=T2                                                             
      T2=T1                                                             
  100 CONTINUE                                                          
      S=S+S-T1                                                          
      BJ=QJ/S                                                           
  300 IF(N) 700,600,600                                                 
  500 DJBES = 0.0D0                                                     
      GO TO 1000                                                        
  600 DJBES = BJ                                                        
      GO TO 1000                                                        
  700 IF( MOD(NN,2) ) 800,600,800                                       
  800 DJBES = -BJ                                                       
      GO TO 1000                                                        
  900 WRITE(6,1001) N,X                                                 
      GO TO 500                                                         
 1000 RETURN                                                            
 1001 FORMAT( 1H ,5X,'THE VALUE OF DJBES IS NOT ACCURATE.  N=',I7,'  , X
     *=',D23.16 )                                                       
      END                                                               
