*** lagged mutual information I(x(t);y(y+tau))
*** input data yy(.,.) in common block yy(1,.) x, yy(2,.) y
*** input parameters:
***** nse: series length (full)
***** kmin,kmax, kstep  range of time lags
***** nq number of equiquantal bins in each variable (the same for both variables)
*** output:
***** r1(.): vector of I(x(t);y(y+tau)), taus are numbered 1,2, ...
***** error: logical variable, if true, problem with binning, results invalid

      subroutine mi2lag(nse,kmin,kmax,kstep,nq,r1,error)

      common /mi2/ yy(2,131072)

      real y1(2,131072)
      real ra(131072),rb(131072),cutp1(1024),r1(256)
      integer mm(1024,1024),m1(1024),m2(1024)
      logical error

      error=.false.

      do iy=1,2

       do jy=1,nse
         ra(jy)=yy(iy,jy)
         rb(jy)=jy
         enddo

      n=nse

      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
          l=l-1
          rra=ra(l)
          rrb=rb(l)
        else
          rra=ra(ir)
          rrb=rb(ir)
          ra(ir)=ra(1)
          rb(ir)=rb(1)
          ir=ir-1
          if(ir.eq.1)then
            ra(1)=rra
            rb(1)=rrb
            goto 33
          endif
        endif
        i=l
        j=l+l
20      if(j.le.ir)then
          if(j.lt.ir)then
            if(ra(j).lt.ra(j+1))j=j+1
          endif
          if(rra.lt.ra(j))then
            ra(i)=ra(j)
            rb(i)=rb(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        go to 20
        endif
        ra(i)=rra
        rb(i)=rrb
      go to 10
 33   continue


*** cutpoints
***** check quanting consistencxy

      xse=nse
      xq=nq
      xx=xse/xq
      nd=xx
        mmx=nd/5
      xx=xse-nd*xq
      if (xx.gt.mmx) then

       write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
       write(*,*)' impossible to quantize:'
       write(*,*)' nse = ',nse
       write(*,*)' nq = ',nq
       write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'

c       stop'error'
       xmi=0.
       error=.true.
       goto 999

       endif


      do j=1,nq-1

       xx=ra(j*nd+1)

       if (ra(j*nd+2).gt.xx) then
        cutp1(j)=xx
        else
        m=0
  110    continue
        m=m+1

        if (m.gt.mmx) then

       write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
       write(*,*)' impossible to equi-quantize:'
       write(*,*)' step ctp1 at'
       write(*,*)' ctp ',j,' of ',nq
       write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'

c       stop'error'
       xmi=0.
       error=.true.
       goto 999

       endif

        xx=ra(j*nd+1+m)

       if (ra(j*nd+2+m).gt.xx) then
        cutp1(j)=xx
        goto 120
        else
        goto 110

        endif

 120     continue
       endif
       enddo

c      write(*,*)' cutp-ed'

***** coding by ctp

        do i=1,nse
         xx=yy(iy,i)

         ix=nq
         do j=1,nq-1

          if ((xx.le.cutp1(j)).and.(ix.eq.nq)) then
           ix=j 

           endif
          enddo
         y1(iy,i)=ix
         enddo

c      write(*,*)' coded'

***** iy cyc:
       enddo
       
       
      nsea=nse-kmax 
      ik=0
      do kk=kmin,kmax,kstep
      ik=ik+1

      do i=1,nq
       m1(i)=0
       m2(i)=0
       do j=1,nq
        mm(i,j)=0
       enddo
       enddo

      pf=float(nsea)

      do i=1,nsea
       i1=y1(1,i)
       i2=y1(2,i+kk)

       mm(i1,i2)=mm(i1,i2)+1
       m1(i1)=m1(i1)+1
       m2(i2)=m2(i2)+1
       enddo

      xmi=0.
      do i=1,nq
       do j=1,nq
        xx=mm(i,j)
        xp=xx/pf
        if (xx.gt.0.) then
         xx=xx/m1(i)
         xx=xx*pf/m2(j)
         xmi=xmi+xp*alog(xx)
        endif
        enddo
       enddo
       
       r1(ik)=xmi
       enddo

 999   continue

      return
      end