c	
c	CROSS.FOR (Cross-correlation of the secondary 
c	event in 1995 and the secondary event in 1996) 
c	QSO 0957+561A,B - April/May 2001 
c
	implicit double precision (a-h),(o-z)
      dimension ta(18),a(18),ea(18),x(18),y(18) 
	dimension tb(18),u(18),v(18),b(18),eb(18)   
	dimension am(18),eam(18),bm(18),ebm(18)
	dimension siga(18),sigb(18),udcf(18,18),tlag(18,18)
	open(unit=1,file='cross15.dat',status='new')
	open(unit=2,file='95s.dat',status='old')
	open(unit=3,file='96s.dat',status='old')
	read(2,*)(ta(i),a(i),ea(i),x(i),y(i),i=1,18)
	read(3,*)(tb(i),u(i),v(i),b(i),eb(i),i=1,18)
	do 10 i=1,18
	if(i.eq.1)am(i)=a(i)
	if(i.gt.1)am(i)=a(i)+am(i-1)
	if(i.eq.1)eam(i)=ea(i)
	if(i.gt.1)eam(i)=ea(i)+eam(i-1)
10	continue
	amedio=am(18)/18.D0
	errora2=(eam(18)/18.D0)**2.D0
	do 20 i=1,18
	if(i.eq.1)bm(i)=b(i)
	if(i.gt.1)bm(i)=b(i)+bm(i-1)
	if(i.eq.1)ebm(i)=eb(i)
	if(i.gt.1)ebm(i)=eb(i)+ebm(i-1)
20	continue
	bmedio=bm(18)/18.D0
	errorb2=(ebm(18)/18.D0)**2.D0
	do 30 i=1,18
	if(i.eq.1)siga(i)=(a(i)-amedio)**2.D0
	if(i.gt.1)siga(i)=(a(i)-amedio)**2.D0+siga(i-1)
30	continue
	sigmaa2=siga(18)/17.D0
	do 40 i=1,18
	if(i.eq.1)sigb(i)=(b(i)-bmedio)**2.D0
	if(i.gt.1)sigb(i)=(b(i)-bmedio)**2.D0+sigb(i-1)
40	continue
	sigmab2=sigb(18)/17.D0
	raiz=dsqrt((sigmaa2-errora2)*(sigmab2-errorb2))
	do 50 i=1,18
	do 60 j=1,18
	udcf(i,j)=(a(i)-amedio)*(b(j)-bmedio)/raiz
	tlag(i,j)=tb(j)-ta(i)
60	continue
50	continue
c	***********************************************************
c	BINS WITH SIZE OF 2*alfa
c	***********************************************************
	alfa=7.5D0
c	***********************************************************
	do 70 k=1,71
	tau=384.D0+k
	bdcf=0.D0
	bdat=0.D0
	do 80 i=1,18
	do 90 j=1,18
	if(tlag(i,j).ge.(tau-alfa).and.tlag(i,j).lt.(tau+alfa))
     *bdcf=bdcf+udcf(i,j)		
      if(tlag(i,j).ge.(tau-alfa).and.tlag(i,j).lt.(tau+alfa))
     *bdat=bdat+1.D0
90	continue
80	continue
	if(bdat.ge.1.D0)write(1,200)tau,bdcf/bdat
70	continue
	stop
200	format(E10.3,4x,E10.3)
	end