C
C     SIMULATION OF N-ROOM QUANTUM MECHANICS
C
C (C) TOSHIAKI IITAKA 1994
C REFERENCE: T.Iitaka, Phys. Rev. E49 (1994) 4684.
C
      PROGRAM NROOM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(PI=3.14159265358979323)
      PARAMETER(NMAX=100)
      COMPLEX*16 P(NMAX),Q(NMAX),R(NMAX),PH0(NMAX),IU,OVR
      DIMENSION HA(NMAX),HB(NMAX-1),V(NMAX)
      DIMENSION HAN(NMAX),HBN(NMAX-1)
C
C     INITIAL DATA
C
      OPEN(20,FILE='tst2.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        V(N)=0
C        HA(N) = V(N) + 1/DX2
      HA(N)=0
   50 CONTINUE
      DO 51 N=1,NMAX-1
          HB(N) = -0.5/DX2
   51 CONTINUE
C
C     SETTING INITIAL WAVE FUNCTION
C
      CALL EIGEN(M,DX,PH0,EM)
      DO 80 N=1,NMAX
      P(N)=PH0(N)
      Q(N)=PH0(N)*EXP(-IU*EM*DT)
   80 CONTINUE
C
C     TIME EVOLUTION LOOP 1000
C
      DO 1000 ITIME=1,MTIME,3
      T= DT*ITIME
C
C      OUTPUT WAVE FUNCTION
C
      IF(MOD(ITIME,3).EQ. 1) THEN
      ANORM=0
      OVR=0.0
      DO 150 N=1,NMAX
        ANORM=ANORM+ABS(Q(N))**2
        OVR=OVR+DCONJG(Q(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/6.* ITIME
      WRITE(20,100) ITIME,EMAX*T,ABS(ANORM-1),
     &ABS(1- ABS(OVR)),ABS(DATAN2(DIMAG(OVR),DBLE(OVR)))
     &,ABS(EM*DT)**3/6. * ITIME 
  100 FORMAT(1H ,I4,30E15.7)
      ENDIF
C
C     TIME EVOLUTION 
C
      R(1) = -2*IU*DT*(HA(1)*Q(1)+HB(1)*Q(2))+P(1)
      DO 70 N=2,NMAX-1
      R(N)=-2*IU*DT*(HB(N-1)*Q(N-1)+HA(N)*Q(N)+HB(N)*Q(N+1))+P(N)
   70 CONTINUE
      R(NMAX)=-2*IU*DT*(HB(NMAX-1)*Q(NMAX-1)+HA(NMAX)*Q(NMAX))+P(NMAX)
C
      P(1) = -2*IU*DT*(HA(1)*R(1)+HB(1)*R(2))+Q(1)
      DO 71 N=2,NMAX-1
      P(N)=-2*IU*DT*(HB(N-1)*R(N-1)+HA(N)*R(N)+HB(N)*R(N+1))+Q(N)
   71 CONTINUE
      P(NMAX)=-2*IU*DT*(HB(NMAX-1)*R(NMAX-1)+HA(NMAX)*R(NMAX))+Q(NMAX)
C
      Q(1) = -2*IU*DT*(HA(1)*P(1)+HB(1)*P(2))+R(1)
      DO 72 N=2,NMAX-1
      Q(N)=-2*IU*DT*(HB(N-1)*P(N-1)+HA(N)*P(N)+HB(N)*P(N+1))+R(N)
   72 CONTINUE
      Q(NMAX)=-2*IU*DT*(HB(NMAX-1)*P(NMAX-1)+HA(NMAX)*P(NMAX))+R(NMAX)
 1000 CONTINUE
      CLOSE(20)
      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