*** (lagged) conditional mutual information for phases with 1d condition
***  I(phi1(t);(phi2(t+tau)-phi2(t))|phi2(t))
*** input data yy(.,.) in common block yy(1,.) phi1, yy(2,.) phi2
*** 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),
***** HERE ONLY 2,4,8 ALLOWED!!
*** output:
***** r3(.,.): CMI as function of tau, taus are numbered 1,2, ...
***** r3(1,.) I(phi2(t);(phi1(t+tau)-phi1(t))|phi1(t)), i.e. 2->1
***** r3(2,.) I(phi1(t);(phi2(t+tau)-phi2(t))|phi2(t)), i.e. 1->2
***** r4(.,.) lagged MI of phases as function of tau, taus are numbered 1,2, ...,
***** but starting from zero lag, i.e. kmin-kstep
****  r4(1,.) I(phi1(t);(phi2(t+tau))
****  r4(2,.) I(phi2(t);(phi1(t+tau))
***** error: logical variable, if true, problem with binning, results invalid


      subroutine c21dfL(nse,kmin,kmax,kstep,nq,r3,r4,error)

      common /cirt/ yy(2,131072)

      real y0(2,131072),yd(2,131072),ctpd(2,16)
      dimension ra(131072),rb(131072),ctp(2,16),
     &r3(2,256),r2(2,256),r4(2,256)

      integer mmm(16,16,16),m12(16,16),m13(16,16),m23(16,16),
     & m1(16),m2(16),m3(16),m4(16),m14(16,16),ival(131072)
       logical error

      error=.false.
      pi=6.*asin(0.5)


      xxq=nq
      h00=alog(xxq)

***** sorting
      do ivar=1,2

      do i=1,nse
        ra(i)=yy(ivar,i)
        rb(i)=i
        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 100

          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
***********************
 100     continue
***********************
***** check quanting consistency

      xse=nse
      xq=nq
      xx=xse/xq
      nd=xx
      xx=nd*xq
      xse1=xse+float(nd)*0.2
      xse2=xse-float(nd)*0.2

      if ((xx.lt.xse2).or.(xx.gt.xse1)) then

       write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
       write(*,*)' impossible to equiquantize:'
       write(*,*)' nse = ',nse,xx
       write(*,*)' nq = ',nq
       write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
       error=.true.
       return
       endif

      do j=1,nq-1

       xx=ra(j*nd+1)

       if (ra(j*nd+2).gt.xx) then
        ctp(ivar,j)=xx
        else
        mm=nd/5
        m=0
  1     continue
        m=m+1

        if (m.gt.mm) then

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

c       stop'error'
       error=.true.
       return

       endif

        xx=ra(j*nd+1+m)
       if (ra(j*nd+2+m).gt.xx) then
        ctp(ivar,j)=xx
        goto 2

        endif

       xx=ra(j*nd+1-m)
       if (ra(j*nd+2-m).gt.xx) then
        ctp(ivar,j)=xx
        goto 2
        else
        goto 1
        endif

 2      continue
       endif
       enddo

**********************************

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

         ix=nq
         do j=1,nq-1

          if ((xx.le.ctp(ivar,j)).and.(ix.eq.nq)) then
           ix=j 

           endif
          enddo
         y0(ivar,i)=ix
         enddo

        enddo  !ivar

**********************************
       ncut=0
       lagov=(kmax-kmin+kstep)/kstep

       nsea=nse-(kmax)-2*ncut

       if (kmin.lt.0) then
       nsea=nse-(kmax+iabs(kmin))-2*ncut
       endif

       pfa=float(nsea)


        i0=1
        if (kmin.lt.0) i0=1-kmin
        i99=0
        if (kmax.gt.0) i99=kmax
**********
********** lag cyc
**********
       ik=0
       do kk=kmin,kmax,kstep
        ik=ik+1
***
*** aby mi faz obsahovala nulovy lag:
        k4=kk
        if (kmin.gt.0) k4=kk-kstep

*********************** diferences, quant
       do ivar=1,2
       do i=1,nsea
        xx=yy(ivar,i+kk)-yy(ivar,i)
        if (xx.lt.-pi) then
c         write(*,*)'dif < -pi ',xx
         xx=xx+2.*pi
c         write(*,*)'                ',xx
         endif
        yd(ivar,i)=xx
        enddo
        enddo

***** sorting
      do ivar=1,2

      do i=1,nsea
        ra(i)=yd(ivar,i)
        rb(i)=i
        enddo

      n=nsea

      l=n/2+1
      ir=n
60    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 150

          endif
        endif
        i=l
        j=l+l
70      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 70
        endif
        ra(i)=rra
        rb(i)=rrb
      go to 60
***********************
 150     continue
***********************
** cutpoints with possible same values
** nq must be power of 2

      xq=nq

      pq=alog(xq)/alog(2.)
      iq=pq
c      if (nq.ne.(2.**iq)) then
c        write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
c        write(*,*)'nq not power of 2 ', nq,2**iq
c        pause
c        endif

*** natural ctp
      nval=0
      do i=2,nsea
       if (ra(i-1).lt.ra(i)) then
        nval=nval+1
        ival(nval)=i
        endif
       enddo
       nval=nval+1

        if (nq.gt.nval) then

       write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
       write(*,*)' impossible to equi-quantize:'
       write(*,*)' nq > # of distinct values'
       write(*,*)  nq,nval
       write(*,*)' DIF faz var, lag = ',ivar,kk
       write(*,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'

       error=.true.
       return

       endif

**  halving
       ipul=nsea/2
       jpul=nq/2  !index half
       ix=1
       xx=ipul-ival(1)
       xx=abs(xx)
       do i=2,nval-1
        xxx=ipul-ival(i)
        xxx=abs(xxx)
        if (xxx.lt.xx) then
         xx=xxx
         ix=i
         endif
        enddo
        ctpd(ivar,jpul)=ix

** quarters
        ict=nsea/4
        jct=nq/4

** l
        ip=ipul-ict
        jp=jpul-jct     
       ix=1
       xx=ip-ival(1)
       xx=abs(xx)
       i9=ctpd(ivar,jpul)-1
       do i=2,i9
        xxx=ip-ival(i)
        xxx=abs(xxx)
        if (xxx.lt.xx) then
         xx=xxx
         ix=i
         endif
        enddo
        ctpd(ivar,jp)=ix

** r
        ip=ipul+ict
        jp=jpul+jct     
       ix=1
       xx=ip-ival(1)
       xx=abs(xx)
       i1=ctpd(ivar,jpul)+1
       do i=i1,nval-1
        xxx=ip-ival(i)
        xxx=abs(xxx)
        if (xxx.lt.xx) then
         xx=xxx
         ix=i
         endif
        enddo
        ctpd(ivar,jp)=ix

*** eighths
      if (nq.ge.8) then
       io=nsea/8
       jo=nq/8
**LL8
        ip=ipul-ict-io
        jp=jpul-jct-jo     
       ix=1
       xx=ip-ival(1)
       xx=abs(xx)
       i9=ctpd(ivar,jpul-jct)-1
       do i=2,i9
        xxx=ip-ival(i)
        xxx=abs(xxx)
        if (xxx.lt.xx) then
         xx=xxx
         ix=i
         endif
        enddo
        ctpd(ivar,jp)=ix

**LP8
        ip=ipul-ict+io
        jp=jpul-jct+jo     
       ix=1
       xx=ip-ival(1)
       xx=abs(xx)
       i1=ctpd(ivar,jpul-jct)+1
       i9=ctpd(ivar,jpul)-1
       do i=i1,i9
        xxx=ip-ival(i)
        xxx=abs(xxx)
        if (xxx.lt.xx) then
         xx=xxx
         ix=i
         endif
        enddo
        ctpd(ivar,jp)=ix

**PL8
        ip=ipul+ict-io
        jp=jpul+jct-jo     
       ix=1
       xx=ip-ival(1)
       xx=abs(xx)
       i1=ctpd(ivar,jpul)+1
       i9=ctpd(ivar,jpul+jct)-1
       do i=i1,i9
        xxx=ip-ival(i)
        xxx=abs(xxx)
        if (xxx.lt.xx) then
         xx=xxx
         ix=i
         endif
        enddo
        ctpd(ivar,jp)=ix

**PP8
        ip=ipul+ict+io
        jp=jpul+jct+jo     
       ix=1
       xx=ip-ival(1)
       xx=abs(xx)
       i1=ctpd(ivar,jpul+jct)+1
       i9=nval-1
       do i=i1,i9
        xxx=ip-ival(i)
        xxx=abs(xxx)
        if (xxx.lt.xx) then
         xx=xxx
         ix=i
         endif
        enddo
        ctpd(ivar,jp)=ix

      endif !nq>=8

*** 16 ??? next time

** prevod indexov na hodnoty
      do i=1,nq
       ix=ctpd(ivar,i)
       jx=ival(ix)
       ctpd(ivar,i)=ra(jx)
       enddo
**********************************

        do i=1,nsea
         xx=yd(ivar,i)

         ix=nq
         do j=1,nq-1

          if ((xx.le.ctpd(ivar,j)).and.(ix.eq.nq)) then
           ix=j 

           endif
          enddo
         yd(ivar,i)=ix
         enddo

        enddo  !ivar

*************** diffs quant

        do i1=1,nq
         do i2=1,nq
           do i3=1,nq
            mmm(i1,i2,i3)=0
            enddo
           m12(i1,i2)=0
           m23(i1,i2)=0
           m13(i1,i2)=0
           m14(i1,i2)=0
           enddo
          m1(i1)=0
          m2(i1)=0
          m3(i1)=0
          m4(i1)=0
          enddo


        do i=ncut+i0,nse-i99-ncut
         i1=y0(1,i)
         i2=yd(1,i)
         i3=y0(2,i)
         i4=y0(2,i+k4)
         

         m1(i1)=m1(i1)+1
         m2(i2)=m2(i2)+1
         m3(i3)=m3(i3)+1
         m4(i4)=m4(i4)+1

         m12(i1,i2)=m12(i1,i2)+1
         m13(i1,i3)=m13(i1,i3)+1
         m23(i2,i3)=m23(i2,i3)+1
         m14(i1,i4)=m14(i1,i4)+1

         mmm(i1,i2,i3)=mmm(i1,i2,i3)+1
         enddo

        xxs=0.
        x1s=0.
        x4s=0.
        do i1=1,nq
           x1=float(m1(i1))/pfa
         do i2=1,nq
           x2=float(m2(i2))/pfa
           x12=float(m12(i1,i2))/pfa
          do i3=1,nq
           x3=float(m3(i3))/pfa

           x13=float(m13(i1,i3))/pfa
           x23=float(m23(i2,i3))/pfa

           xxx=float(mmm(i1,i2,i3))/pfa

           if (xxx.gt.0.) then

             xxs=xxs+xxx*alog((xxx*x1)/(x12*x13))
             endif

            enddo


             x14=float(m14(i1,i2))/pfa
             x4=float(m4(i2))/pfa
           if (x14.gt.0.) then
             x4s=x4s+x14*alog(x14/(x1*x4))
             endif


           enddo
          enddo

         r3(1,ik)=xxs

         r4(1,ik)=x4s
*** opacne:

        do i1=1,nq
         do i2=1,nq
           do i3=1,nq
            mmm(i1,i2,i3)=0
            enddo
           m12(i1,i2)=0
           m23(i1,i2)=0
           m13(i1,i2)=0
           m14(i1,i2)=0
           enddo
          m1(i1)=0
          m2(i1)=0
          m3(i1)=0
          m4(i1)=0
          enddo


        do i=ncut+i0,nse-i99-ncut
         i1=y0(2,i)
         i2=yd(2,i)
         i3=y0(1,i)
         i4=y0(1,i+k4)
         

         m1(i1)=m1(i1)+1
         m2(i2)=m2(i2)+1
         m3(i3)=m3(i3)+1
         m4(i4)=m4(i4)+1

         m12(i1,i2)=m12(i1,i2)+1
         m13(i1,i3)=m13(i1,i3)+1
         m23(i2,i3)=m23(i2,i3)+1
         m14(i1,i4)=m14(i1,i4)+1

         mmm(i1,i2,i3)=mmm(i1,i2,i3)+1
         enddo


        xxs=0.
        x1s=0.
        x4s=0.
        do i1=1,nq
           x1=float(m1(i1))/pfa
         do i2=1,nq
           x2=float(m2(i2))/pfa
           x12=float(m12(i1,i2))/pfa
          do i3=1,nq
           x3=float(m3(i3))/pfa

           x13=float(m13(i1,i3))/pfa
           x23=float(m23(i2,i3))/pfa

           xxx=float(mmm(i1,i2,i3))/pfa

           if (xxx.gt.0.) then

             xxs=xxs+xxx*alog((xxx*x1)/(x12*x13))
             endif

            enddo


             x14=float(m14(i1,i2))/pfa
             x4=float(m4(i2))/pfa
           if (x14.gt.0.) then
             x4s=x4s+x14*alog(x14/(x1*x4))
             endif


           enddo
          enddo

         r3(2,ik)=xxs

         r4(2,ik)=x4s

         enddo !lag cyc


      return
      end