PROGRAM SYNCHRO
*	Given a bivariate input data (infile) of NDATOS data points, the program 
*	calculates the cross-correlation and the non-linear interpendences S,H and N.
*	Results are ploted in the screen with the following order:
*	S(X|Y) S(Y|X) H(X|Y) H(Y|X) N(X|Y) N(Y|X) cross-correlation.
*	Parameters to be set are: m (embedding dimension), tau (time delay), 
*	theiler (theiler correction in # data points) and NN (# of nearest neighbors).
*

      parameter (NDATOS=1000)
      parameter ( m=10, tau=2, theiler=50, NN = 10)
      REAL * 8 X1(NDATOS), x1n(NDATOS), Y1(NDATOS), y1n(NDATOS)
      REAL * 8 AUXX(1001), AUXY(1001), INDEXX(1001), INDEXY(1001)
      REAL * 8 DISTX(NDATOS), DISTY(NDATOS),HXY,HYX
      REAL * 8  xmean, ymean, xvar, yvar, aa, bb,NXY,NYX
      REAL * 8 RXX, RYY, RXY, RYX, SXY, SYX, RRX, RRY, cross
      character * 30 infile
*      
      infile = 'example_a.asc'
      open(unit = 200, file = infile, status = 'old')
             xmean = 0
             ymean = 0
         do i = 1, ndatos                                   ! read data
             read (200,*) aa,bb
             x1(i) = aa
             y1(i) = bb
             xmean = xmean + x1(i)
             ymean = ymean + y1(i)
         end do
             xmean = xmean / ndatos
             ymean = ymean / ndatos
             xvar = 0
             yvar = 0
           do i = 1, ndatos                          !normalization (x1n ; y1n)
             xvar = xvar + ( x1(i) - xmean )**2
             yvar = yvar + ( y1(i) - ymean )**2
           end do
             xvar = sqrt (xvar / (ndatos-1))
             yvar = sqrt (yvar / (ndatos-1))
             cross = 0
           do i= 1, ndatos
             x1n(i) = (x1(i) - xmean) / xvar
             y1n(i) = (y1(i) - ymean) / yvar
             cross = cross + x1n(i)*y1n(i)                    !cross correlation
           end do  
             cross = cross/ndatos
        close(unit = 200)  
*
	SXY = 0
	SYX = 0
	HXY = 0
	HYX = 0
	NXY = 0
	NYX = 0
      DO I = 1, NDATOS-(m-1)*tau-1
                do k=1,nn                                      !INICIALIZE AUX
                   AUXX(K) = 100000000
                    INDEXX(K) = 100000000
                    AUXY(K) = 100000000
                   INDEXY(K) = 100000000
                end do
                   AUXX(NN+1) = 0.D0
                   AUXY(NN+1) = 0.D0
                   INDEXX(NN+1) = 100000000
                   INDEXY(NN+1) = 100000000
                   RRX = 0
                   RRY = 0
         DO J = 1, NDATOS-(m-1)*tau-1                         !*****************
               DISTX(J) = 0
               DISTY(J) = 0
          do k=0,m-1                                            !DISTANCES
           DISTX(J) = DISTX(J)+(X1(I+k*tau)-X1(J+k*tau))**2
           DISTY(J) = DISTY(J)+(Y1(I+k*tau)-Y1(J+k*tau))**2
          end do
*
           if ((abs(i-j)) .gt. theiler) then  
              IF (DISTX(J) .LT. AUXX(1)) THEN
                 DO K=1,NN+1
                    IF (DISTX(J) .LT. AUXX(K)) THEN
                       AUXX(K) = AUXX(K+1)
                       INDEXX(K) = INDEXX(K+1)
                      ELSE
                       AUXX(K-1) = DISTX(J)
                       INDEXX(K-1) = J
                       goto 10
                    END IF
                 END DO
              END IF
 10           continue
              IF (DISTY(J) .LT. AUXY(1)) THEN
                 DO K=1,NN+1
                    IF (DISTY(J) .LT. AUXY(K)) THEN
                       AUXY(K) = AUXY(K+1)
                       INDEXY(K) = INDEXY(K+1)
                      ELSE
                       AUXY(K-1) = DISTY(J)
                       INDEXY(K-1) = J
                       goto 20
                    END IF
                 END DO
              END IF          
 20           continue
           end if    
              RRX = RRX + DISTX(J)        !SIZE OF THE ATTRACTORS
              RRY = RRY + DISTY(J)
         END DO                                               !***************** 
              RXX = 0
              RYY = 0
              RXY = 0
              RYX = 0
           DO K=1,NN
             RXX = AUXX(K) + RXX
              RYY = AUXY(K) + RYY
              RXY = DISTX(INDEXY(K)) + RXY
             RYX = DISTY(INDEXX(K)) + RYX
           END DO
              RXX = RXX/NN
              RYY = RYY/NN
              RXY = RXY/NN
              RYX = RYX/NN
             SXY = SXY + RXX/RXY
             SYX = SYX + RYY/RYX
             HXY =  HXY + LOG(RRX/((NDATOS-(M-1)*tau-2)) / RXY)
             HYX =  HYX + LOG(RRY/((NDATOS-(M-1)*tau-2)) / RYX) 
             NXY =  NXY + RXY / (RRX/((NDATOS-(M-1)*tau-2))) 
             NYX =  NYX + RYX / (RRY/((NDATOS-(M-1)*tau-2)))  
      END DO
             SXY = SXY/(NDATOS-(M-1)*tau-1)
             SYX = SYX/(NDATOS-(M-1)*tau-1)
             HXY = HXY/(NDATOS-(M-1)*tau-1)
             HYX = HYX/(NDATOS-(M-1)*tau-1)
             NXY = 1 - NXY/(NDATOS-(M-1)*tau-1)
             NYX = 1 - NYX/(NDATOS-(M-1)*tau-1)
*              
             write(6,61) SXY, SYX, HXY, HYX, NXY, NYX, cross
61     format (7(x,e12.5)) 
      END