ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c ERRORS5SIM.FOR (Error in the measurement of the delay c c through the secondary 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(18),ao(18),ea(18),x(18),y(18) dimension tb(18),u(18),v(18),bo(18),eb(18) dimension a(18),b(18) dimension am(18),bm(18),siga(18),sigb(18) dimension udcf(18,18),tlag(18,18),udcfa(18,18),tlaga(18,18) dimension tauc(71),taua(71),cc(71),ac(71) dimension teta(41),test(41) open(unit=1,file='errors5sim.dat',status='new') open(unit=2,file='95s.dat',status='old') open(unit=3,file='96s.dat',status='old') read(2,*)(ta(i),ao(i),ea(i),x(i),y(i),i=1,18) read(3,*)(tb(i),u(i),v(i),bo(i),eb(i),i=1,18) call g05cbf(0) c *********************************************************** c BINS WITH SIZE OF 2*alfa c *********************************************************** alfa=2.5D0 c cont425=0.D0 cont426=0.D0 cont427=0.D0 cont428=0.D0 cont429=0.D0 cont430=0.D0 cont431=0.D0 cont432=0.D0 cont433=0.D0 cont434=0.D0 cont435=0.D0 cont436=0.D0 cont437=0.D0 cont438=0.D0 cont439=0.D0 cont440=0.D0 do 3 n=1,100000 c *********************************************************** c SYNTHETIC LIGHT CURVES AND THEIR PROPERTIES c *********************************************************** do 5 i=1,18 a(i)=g05ddf(ao(i),ea(i)) 5 continue do 7 i=1,18 b(i)=g05ddf(bo(i),eb(i)) 7 continue do 10 i=1,18 if(i.eq.1)am(i)=a(i) if(i.gt.1)am(i)=a(i)+am(i-1) 10 continue amedio=am(18)/18.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) 20 continue bmedio=bm(18)/18.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*sigmab2) raiza=sigmaa2 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 do 70 i=1,18 do 80 j=1,18 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,71 tau=384.D0+k bdcf=0.D0 bdat=0.D0 do 100 i=1,18 do 110 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 110 continue 100 continue tauc(k)=tau cc(k)=bdcf/bdat 90 continue c *********************************************************** c AUTOCORRELATION c *********************************************************** do 120 k=1,71 tau=-36.D0+k bdcf=0.D0 bdat=0.D0 do 130 i=1,18 do 140 j=1,18 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,41 teta(k)=399.D0+k delta=0.D0 adat=0.D0 do 160 i=1,71 do 170 j=1,71 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,41 testmin=dmin1(test(i),testmin) 180 continue do 190 i=1,41 if(test(i).eq.testmin)delay=teta(i) 190 continue if(delay.eq.425.D0)cont425=cont425+1.D0 if(delay.eq.426.D0)cont426=cont426+1.D0 if(delay.eq.427.D0)cont427=cont427+1.D0 if(delay.eq.428.D0)cont428=cont428+1.D0 if(delay.eq.429.D0)cont429=cont429+1.D0 if(delay.eq.430.D0)cont430=cont430+1.D0 if(delay.eq.431.D0)cont431=cont431+1.D0 if(delay.eq.432.D0)cont432=cont432+1.D0 if(delay.eq.433.D0)cont433=cont433+1.D0 if(delay.eq.434.D0)cont434=cont434+1.D0 if(delay.eq.435.D0)cont435=cont435+1.D0 if(delay.eq.436.D0)cont436=cont436+1.D0 if(delay.eq.437.D0)cont437=cont437+1.D0 if(delay.eq.438.D0)cont438=cont438+1.D0 if(delay.eq.439.D0)cont439=cont439+1.D0 if(delay.eq.440.D0)cont440=cont440+1.D0 3 continue write(1,*)'100000 SIMULATIONS (second-second)' write(1,*)'*************************' write(1,*)'#delays-425 =',cont425 write(1,*)'#delays-426 =',cont426 write(1,*)'#delays-427 =',cont427 write(1,*)'#delays-428 =',cont428 write(1,*)'#delays-429 =',cont429 write(1,*)'#delays-430 =',cont430 write(1,*)'#delays-431 =',cont431 write(1,*)'#delays-432 =',cont432 write(1,*)'#delays-433 =',cont433 write(1,*)'#delays-434 =',cont434 write(1,*)'#delays-435 =',cont435 write(1,*)'#delays-436 =',cont436 write(1,*)'#delays-437 =',cont437 write(1,*)'#delays-438 =',cont438 write(1,*)'#delays-439 =',cont439 write(1,*)'#delays-440 =',cont440 write(1,*)'#delays-(425,440) =',cont425+cont426+cont427+ *cont428+cont429+cont430+cont431+cont432+cont433+ *cont434+cont435+cont436+cont437+cont438+cont439+cont440 stop end