program cross

*     Program to generate cross correlation with and with out
*     noise
      integer*4 mxn
      parameter ( mxn = 4096)   ! Maximum number of values in TS

      real*4 x(mxn), y(mxn), nx(mxn), ny(mxn), rxx(mxn), rxy(mxn)
      real*4 ss, sn, c, snr, xs
      integer*4 i,j,k, n, ns,ne, num
      character*1 tb      ! tab character of tab-delimited file.
      logical clipping

*     Set the sample size (n) and where we will start correlation (ns)
      n = 1024
      ns = 256
C     n = 4096
C     ns = 256
*     Set sample number of end of correlation (symmetric) and compute
*     the number of samples
      ne = n-ns
      num = ne - ns 

*     Set clipping true to make samples be +-1, set false to keep original
*     samples
      clipping = .true.

*     Set the tab character
      tb = char(9)

*     Set the SNR desired and compute the sigmas of the signal and
*     noise to acheive this SNR (note variance of uniform distribution
*     between -0.5 and 0.5 is 1/12.      
      SNR = 0.1
      if( SNR.lt.1000 ) then
          c = 12.d0/(1+SNR)
      else      ! Take to be infinite SNR
          c = 0.d0
      endif
* Example values
c     c = 10.90909091d0    ! SNR = 0.1
c     c = 6.d0             ! SNR = 1

*     Compute standard deviation of signal (ss) and of noise (sn)
*     (Values choosen so that sum of signals will have variance of 1.0)
      ss = sqrt(12.d0-c)
      sn = sqrt(c)

*     Generate x and y time series.  Y is simple x lagged by 50 steps.
*     Each time series has independent noise.
*     Data is written to unit 20, the auto and cross correlations are 
*     written to unit 21.
      if( clipping ) then
         write(20,'(5a)') 'Sample',tb,'X Clipped',tb,'Y Clipped'
      else
         write(20,'(5a)') 'Sample',tb,'X',tb,'Y'
      endif

*     Loop over n-samples generating the time series
      do i = 1, n
*        Value of signal at this time
         xs = (rand()-0.5d0)*ss 
*        Add the noise to the signal
         x(i) = xs + (rand()-0.5d0)*sn
*        generate y as a lagged version of x with independent noise
         if( i+50.le.n ) then
            y(i+50) = xs + (rand()-0.5d0)*sn
         endif
*        Generate random values for y before the start of the lagged 
*        portion (x also has 50 samples at the end that do not appear
*        in the y-timeseries.
         if( i.lt.50) then
            y(i) = (rand()-0.5d0)*ss + (rand()-0.5d0)*sn
         endif

*        If we are clipping replace values with +-1 depending on sign
         if( clipping ) then
            if ( x(i).ge.0 ) then
                x(i) = 1.0
            else
                x(i) = -1.0
            end if
            if ( y(i).ge.0 ) then
                y(i) = 1.0
            else
                y(i) = -1.0
            end if
         endif
*        Write out the data for plotting (offset values by +-2 so that
*        separate time series can be seen).         
         write(20,100) i, tb, x(i)+2.0,tb,  y(i)-2.0
 100     format(i5,a,F10.6,a,f10.6)
      end do

****  Now form the correlation function
*     Write out the header line for the output to unit 21
      if( sn.gt.0 ) then
          if( clipping ) then
             write(21,'(5a,f6.2)') 
     .          'Lag',tb,'\\fr\\nxx Clip',tb,
     .          '\\fr\\nxy Clip SNR ',snr
          else
             write(21,'(5a,f6.2)') 
     .          'Lag',tb,'\\fr\\nxx',tb,'\\fr\\nxy SNR ',snr
          endif
      else
          if( clipping ) then
             write(21,'(5a)') 
     .          'Lag',tb,'\\fr\\nxx Clip',tb,
     .          '\\fr\\nxy Clip SNR Infinite'
          else
             write(21,'(5a)') 
     .          'Lag',tb,'\\fr\\nxx',tb,'\\fr\\nxy SNR Infinite'
          endif
      end if

*     Now do the correlation just over +-100 laggs about the center
*     (Offset is 50 so will appear half way to one side).
      do k = -100,100
         rxx(k+101) = 0.d0
         rxy(k+101) = 0.d0
         do j = ns,ne
*           Form autocorrelation of x and cross correlation of x and y
            rxx(k+101) = rxx(k+101)+x(j)*x(j+k)
            rxy(k+101) = rxy(k+101)+x(j)*y(j+k)
         end do
*        Nomralize by number of samples
         rxx(k+101) = rxx(k+101)/num
         rxy(k+101) = rxy(k+101)/num
*        Write out the values`
         write(21,210) k,tb,rxx(k+101),tb,rxy(k+101)
 210     format(i5,a,F10.6,a,f10.6)
      end do

      end