      program exact_2N_1D

c     exact transfer matrix product calculation
c     for 2N system in 1D
c     with initial state equal to 1 
c     for all spatial lattice sites

      implicit integer(i-n)
      implicit double precision(a-h,o-y)
      implicit complex*16(z)

ccccccccc

      parameter(Lt = 50)
      parameter(L = 6)
      parameter(ainv = 100.D0)
      parameter(atinv = 100.D0)
      parameter(alphat = ainv/atinv)
      parameter(c = -0.2D0)
      parameter(am = 938.92D0)
      dimension vrel(0:L-1,0:Lt)
      dimension temp(0:L-1)
      dimension overlap(0:Lt)

cccccc

      do nx = 0,L-1
         vrel(nx,0) = 1.D0/dsqrt(1.D0*L)
      enddo     
      overlap(0) = 1.D0

      do nt = 1,Lt
         do nx = 0,L-1
            temp(nx) = 
     $           vrel(nx,nt-1)*(1.D0 - 1.D0/am*alphat)
     $           + vrel(mod(nx-1+L,L),nt-1)/(2.D0*am)*alphat
     $           + vrel(mod(nx+1,L),nt-1)/(2.D0*am)*alphat
         enddo
         do nx = 0,L-1
            vrel(nx,nt) = 
     $           temp(nx)*(1.D0 - 1.D0/am*alphat)
     $           + temp(mod(nx-1+L,L))/(2.D0*am)*alphat
     $           + temp(mod(nx+1,L))/(2.D0*am)*alphat
         enddo
         vrel(0,nt) = vrel(0,nt) - c*alphat*vrel(0,nt-1)
         overlap(nt) = 0.D0
         do nx = 0,L-1
            overlap(nt) = overlap(nt) + vrel(nx,nt)/dsqrt(1.D0*L)
         enddo         
      enddo

      do nt = 1,Lt
         ratio = overlap(nt)/overlap(nt-1)
         write(*,*)'nt = ',nt,
     $        'energy = ', -dlog(ratio)*atinv,'(MeV)'
      enddo

      stop
      end

