ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc	
c                                                                 c
c	ERRORP5SIM.FOR (Error in the measurement of the delay       c
c	through the main events. We use bins with size of 5         c
c	days and do simulations of pairs of possible true           c
c	signals) - NAG Fortran Library                              c
c	February 2002                                               c
c	                                                            c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	implicit double precision (a-h),(o-z)
      dimension ta(13),ao(13),ea(13),x(13),y(13) 
	dimension tb(14),u(14),v(14),bo(14),eb(14)
	dimension a(13),b(14)   
	dimension am(13),bm(14),siga(13),sigb(14)
	dimension udcf(13,14),tlag(13,14),udcfa(13,13),tlaga(13,13)
	dimension tauc(201),taua(201),cc(201),ac(201)
	dimension teta(161),test(161)
	open(unit=1,file='errorp5sim.dat',status='new')
	open(unit=2,file='95p.dat',status='old')
	open(unit=3,file='96p.dat',status='old')
	read(2,*)(ta(i),ao(i),ea(i),x(i),y(i),i=1,13)
	read(3,*)(tb(i),u(i),v(i),bo(i),eb(i),i=1,14)
	call g05cbf(0)
c	***********************************************************
c	BINS WITH SIZE OF 2*alfa 
c	***********************************************************
	alfa=2.5D0
c
	cont416=0.D0
	cont41625=0.D0
	cont41650=0.D0
	cont41675=0.D0
	cont417=0.D0
	cont41725=0.D0
	cont41750=0.D0
	cont41775=0.D0
	cont418=0.D0
	do 3 n=1,100000
c	***********************************************************
c	SYNTHETIC LIGHT CURVES AND THEIR PROPERTIES
c	***********************************************************
	do 5 i=1,13
	a(i)=g05ddf(ao(i),ea(i))
5	continue
	do 7 i=1,14
	b(i)=g05ddf(bo(i),eb(i))
7	continue
	do 10 i=1,13
	if(i.eq.1)am(i)=a(i)
	if(i.gt.1)am(i)=a(i)+am(i-1)
10	continue
	amedio=am(13)/13.D0
	do 20 i=1,14
	if(i.eq.1)bm(i)=b(i)
	if(i.gt.1)bm(i)=b(i)+bm(i-1)
20	continue
	bmedio=bm(14)/14.D0
	do 30 i=1,13
	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(13)/12.D0
	do 40 i=1,14
	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(14)/13.D0
	raiz=dsqrt(sigmaa2*sigmab2)
	raiza=sigmaa2
	do 50 i=1,13
	do 60 j=1,14
	udcf(i,j)=(a(i)-amedio)*(b(j)-bmedio)/raiz
	tlag(i,j)=tb(j)-ta(i)
60	continue
50	continue
	do 70 i=1,13
	do 80 j=1,13
	udcfa(i,j)=(a(i)-amedio)*(a(j)-amedio)/raiza
	tlaga(i,j)=ta(j)-ta(i)
80	continue
70	continue
c	***********************************************************
c	CROSS-CORRELATION
c	***********************************************************
	do 90 k=1,201
	tau=394.75D0+k*0.25D0
	bdcf=0.D0
	bdat=0.D0
	do 100 i=1,13
	do 110 j=1,14
	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
110	continue
100	continue
	tauc(k)=tau
	cc(k)=bdcf/bdat
90	continue
c	***********************************************************
c	AUTOCORRELATION 
c	***********************************************************
	do 120 k=1,201
	tau=-25.25D0+k*0.25D0
	bdcf=0.D0
	bdat=0.D0
	do 130 i=1,13
	do 140 j=1,13
	if(tlaga(i,j).ge.(tau-alfa).and.tlaga(i,j).lt.(tau+alfa))
     *bdcf=bdcf+udcfa(i,j)		
      if(tlaga(i,j).ge.(tau-alfa).and.tlaga(i,j).lt.(tau+alfa))
     *bdat=bdat+1.D0
140	continue
130	continue
	taua(k)=tau
	ac(k)=bdcf/bdat
120	continue
c	***********************************************************
c	delta^2 TEST  
c	***********************************************************
	do 150 k=1,161
	teta(k)=399.75D0+k*0.25D0
	delta=0.D0
	adat=0.D0
	do 160 i=1,201
	do 170 j=1,201
	if(taua(j).eq.(tauc(i)-teta(k)))delta=delta+
     *(cc(i)-ac(j))**2.D0
	if(taua(j).eq.(tauc(i)-teta(k)))adat=adat+1.D0
170	continue
160	continue
	test(k)=delta/adat 
150	continue
	testmin=test(1)
	do 180 i=2,161
	testmin=dmin1(test(i),testmin)
180 	continue
	do 190 i=1,161
	if(test(i).eq.testmin)delay=teta(i)
190	continue
	if(delay.eq.416.D0)cont416=cont416+1.D0
	if(delay.eq.416.25D0)cont41625=cont41625+1.D0
	if(delay.eq.416.50D0)cont41650=cont41650+1.D0
	if(delay.eq.416.75D0)cont41675=cont41675+1.D0
	if(delay.eq.417.D0)cont417=cont417+1.D0
	if(delay.eq.417.25D0)cont41725=cont41725+1.D0
	if(delay.eq.417.50D0)cont41750=cont41750+1.D0
	if(delay.eq.417.75D0)cont41775=cont41775+1.D0
	if(delay.eq.418.D0)cont418=cont418+1.D0
3	continue
	write(1,*)'100000 SIMULATIONS (main-main)'
	write(1,*)'**********************************'
	write(1,*)'#delays-416 =',cont416
	write(1,*)'#delays-41625 =',cont41625
	write(1,*)'#delays-41650 =',cont41650
	write(1,*)'#delays-41675 =',cont41675
	write(1,*)'#delays-417 =',cont417
	write(1,*)'#delays-41725 =',cont41725
	write(1,*)'#delays-41750 =',cont41750
	write(1,*)'#delays-41775 =',cont41775
	write(1,*)'#delays-418 =',cont418
	stop
	end